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THE SELF-ADAPTIVE GRID CODE, SAGEv2 


Carol B. Davies* and Ethiraj Venkatapathyt 
Ames Research Center 

SUMMARY 

The original SAGE code (Version 1) is described in NASA TM 103905 (1992). This new 
report on Version 2 includes all the information in the original publication plus all upgrades and 
changes to the SAGE code since that time. The two most significant upgrades are the inclusion of 
a finite-volume option and the ability to adapt and manipulate zonal-matching multiple-grid 
files. In addition, the original SAGE code has been upgraded to Version 1.1 and includes all 
options mentioned in this report, with the exception of the multiple grid option and its 
associated features. Since Version 2 is a larger and more complex code, it is suggested (but not 
required) that Version 1.1 be used for single-grid applications. This document contains all the 
information required to run both versions of SAGE. 

The formulation of the adaption method is described in the first section of this document. 
The second section is presented in the form of a user guide that explains the input and execution 
of the code. The third section provides many examples. 

Successful application of the SAGE code in both two and three dimensions for the solution 
of various flow problems has proven the code to be robust, portable, and simple to use. Although 
the basic formulation follows the method of Nakahashi and Deiwert, many modifications have 
been made to facilitate the use of the self-adaptive grid method for complex grid structures. 
Modifications to the method and the simple but extensive input options make this a flexible and 
user-friendly code. The SAGE code can accommodate two-dimensional and three-dimensional, 
finite-difference and finite-volume, single grid and zonal-matching multiple grid flow problems. 

1. METHOD 
1.1 Introduction 

Solution-adaptive grid methods have become useful tools for efficient and accurate flow 
predictions. In supersonic and hypersonic flows, strong gradient regions such as shocks, contact 
discontinuities, and shear layers require careful distribution of grid points to minimize grid error 
and thus produce accurate flow-field predictions. Frequently, the generation of the first grid 
topology does not adequately capture these flow structures. It has been shown that an effective 
way of obtaining accurate solutions is by intelligently redistributing (i.e., adapting) the original 
grid points based on the initial flow-field solution and then computing a new solution using the 
adapted grid. The cost efficiency of grid adaptions in terms of CPU time depends on the basic 
formulation of the adaptive-grid solver. A short historical review and a list of references on the 
adaptive grid procedure used in the SAGE code is given in Davies and Venkatapathy (1991). 
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The self-adaptive grid procedure outlined by Nakahashi and Deiwert (1985) has evolved 
into a flexible tool for adapting and restructuring both two-dimensional (2-D) and three- 
dimensional (3-D) grids. The adaptive-grid method is based on grid-point redistribution through 
local error minimization. The procedure is analogous to applying tension and torsion spring 
forces proportional to the local flow gradient at every grid point and finding the equilibrium 
position of the resulting system of grid points. The multidimensional problem of grid adaption is 
split into a series of one-dimensional (1-D) problems along the computational coordinate lines. 
The reduced 1-D problem then requires a tridiagonal solver to find the location of grid points 
along a given coordinate line. Multidirectional adaption is achieved by the sequential application 
of the method in each coordinate direction. This approach produces an extremely CPU-efficient 

algorithm. • , 

The tension forces direct the redistribution of points to the strong gradient regions. The 
torsion forces relate information between the families of lines adjacent to one another, to 
maintain smoothness and a measure of orthogonality of grid lines. These smoothness and 
orthogonality constraints are direction-dependent, since they relate only the coordinate lines that 
are being adapted to the neighboring lines that have already been adapted. This implies that the 
solutions are nonunique and depend on the order and direction of adaption. 

The Self-Adaptive Grid codE (SAGE) code has been built with many flexible elements that 
make it user-friendly. The second section of this report is a user guide, which is independent of 
the first section and can be used as a separate document. However, the nomenclature and details 
of the grid adaption procedure within the code consistently follow the analysis, allowing the user 
to understand the code and implement any individual changes, if desired. The user guide 
describes the input parameters and their effects as well as the adaption procedure, execution 
commands, and subroutines, and provides many examples. 

1.2 Formulation of Adaptive-Grid Scheme 

As stated in the introduction, adaption takes place as a series of one-dimensional 
adaptions. Figure 1 illustrates this concept: three constant k planes of an initial grid are shown in 
Fig. 1(a) and a flow-field solution has been obtained using this unadapted grid. The points in this 
grid are then adapted to the flow solution, starting on the first line (j=l) on the lower plane k=l. 
In Fig. 1(b), the first plane has already been adapted and the second plane is the current adaption 
plane. The current adaption line (j) is shown, with previous lines already adapted and 

subsequent lines awaiting adaption. 
The third plane is still in its original 
form. Adaption is performed in this 
line-by-line, plane-by-plane manner 
until all requested planes are 
complete. It is then possible to 
perform an adaption in a second 
direction, "adapting" on top of the 
already adapted grid. The number 
and order of adaptions are arbitrary 
and depend on the type of flow 
problem and the purpose of the 
adaption. 



Figure 1. 3-D adaption ; (a) Initial grid; 
(b) directional adaption. 
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Figure 2 shows a segment 
of the current adaption line in 
more detail. The lower plane 
has already been adapted and the 
upper plane is currently being 
adapted. Four forces control the 
redistribution of points along 
each line: the two tension 
springs that act on each side of a 
node, and the two torsion 
springs that control the 
smoothness of the grid. The 

tension forces co cluster the 
redistributed points into the 
high-gradient regions. The 
torsion forces (x and y) restrain 
the tension forces to maintain 
continuity between sequentially 
adapted lines. As shown in Fig. 

2, the T force acts from the previously adapted line within the current plane and the V|/ force acts 
from the previously adapted line in the preceding plane. 

A 3-D grid is described in terms of its computational coordinate directions (i,j,k) and its 
physical coordinates (x,y,z). An adaption can take place in any or all of the computational 
coordinate directions and the SAGE code permits any combination. However, for clarity, this 
analysis assumes that all adaptions are performed in the i direction and that stepping (within a 
plane surface) occurs in the j direction. Thus k is the marching direction for plane stepping. The 
current adaption line is thus j in the plane k where lines etc., have already been adapted. 

In addition, this report uses the term "plane" to mean any 3-D coordinate surface. 

The first step in the formulation of the adaption algorithm is to consider the adaption of a 
single line, with no torsional constraints 
(i.e., no smoothness control). Figure 3 shows 
a 2-D surface, where the arc length at A (i.e., 
s i jk , with s, =0) along the current adaption 
line, j, is defined as 

s i - *,-i +(y, — y,_, ; 2 +tz j -z,_, J 2 (1) 

A tension force, f (<£>,), is defined to 
act between each i and i+1 node such that 

CO.Av, = K (2) 

where to, the tension constant, is a 
weighting function based on the flow 
gradient. Ay, = y, +1 - s t , and K is the resultant 
force. To redistribute the points along a line 
with the minimum solution error, the 
adaptive procedure defines the weighting 
function as proportional to the derivative of 





Figure 3. Current adaption line j on a 2-D surface. 



to =tension spring between each node on current adaption line 
T =torsion spring within current plane at previously adapted line 
V =torsion spring between planes at current line on already adapted plane 

Figure 2. Line-by-line adaption 
showing tension and torsion forces. 
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the flow variable (Nakahashi and Deiwert 1985). In this formulation, (0 is defined as a function of 
the normalized flow gradient, /, such that 

co, =1 + 4 /* ( 3 ) 

where A and B are constants directly related to the grid spacing and are chosen to maintain the 
grid intervals to within the limits (A s MIN and Aj wax ) set by the user on input. A and B are defined 
in appendix I of this section. 

Equation (2) is written for each node on the line, giving a 1-D formulation that can be 
solved directly for As r Taking the sum of both sides of Eq. (2) gives 

(=i <=i 


giving 


Substituting back in Eq. (2), we obtain 


K = s. 


As, = s m 




f CO 


*“ ©I 


(4) 


(5) 


In the SAGE code, this 1-D solution technique is used along the initial adaption line of the initial 
adaption plane, where no directional information is available. Continuing this approach for 
successive line-by-line adaptions will not create a mesh that is sufficiently smooth for input into 
computational flow-field codes. To ensure a more reasonable grid, the redistribution of points 
(driven by the tension springs) is constrained by torsion terms defined on two adjacent adapted 
lines: one on the current plane and one on the preceding plane. Within the current plane, the 

torsion parameter T represents the magnitude of the torsion force that maintains smoothness and 
orthogonality between the node (i,j,k) and the nodes (i,j-l,k) and (i,j-2,k). This is the only torsion 
parameter available on the initial surface. For subsequent surfaces, the torsion parameter \\i 
constrains the movement between the same node (i,j,k) and the nodes (i,j,k-l) and (i,j,k-2) on 
adjacent computational planes. The torsion terms are evaluated as t,( s , ~ ) and YjfSf 
where the torsion parameters x and y define the magnitude of the torsion forces and 5 ' and s' 
define their inclination (i.e., orthogonality and smoothness). 

To introduce the torsion forces to the system of equations, we first rewrite Eq. (2) to 
represent the force balance, 

to , ( s M -S:)- 00 , ( 5 , - s,_i ) = 0 


and then add the torsion terms to obtain 

which is rearranged to produce the coefficients of the tridiagonal matrix used to solve for s jr 

co,_,s,_, - (to, + 00,-1 + T, + V, )s, + (ti,s l+l = -x/ - s' (8) 

This equation is written for each interval along the adaption line, producing a system of (n l — 1) 
equations. Since 5 , and s„ are known, there are n, - 2 rows in the matrix. This equation is solved 

iteratively (updating co, at new s,) until 10 " 3 This s y stem of e q uations is 
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central to the adaption technique used in the SAGE code and most of the description that follows 
pertains to deriving the coefficients of this equation. 


1.3 Auxiliary Equations 


1.3.1 Tension parameter, co,. As described in the formulation, co, acts as a tension force in 
the interval ( s, +1 - s t ) and can be imagined to be a spring (aligned with the grid line) connecting 
the two grid points, as shown in Fig. 3. This tension parameter is defined in Eq. (2) as 

co, A j, = K 


where /, is a function of the gradient of the flow variable, q: 


f,= 




mm 


f max f mm 


and /,= 


d<7. <7,+|-g, 

ds As j 


(9) 


The standard format of the user's input flow-field file contains the conservative flow 
variables Q, as defined in the plotting software package PLOT3D (1990), which (for a 3-D dataset) 

assumes these flow variables are p, pu, pv, pw, and e. Since the adaption is based on a scalar 
function q, this function can be evaluated as a user-specified combination of flow variables Q. 
Pressure, Mach number and temperature ratio are also included as adaption variables and are 
internally computed, assuming perfect gas relations. However, Q need not be restricted to 
conservative flow variables, but can be any combination of user-specified flow variables that 
represent the flow field. 

To remove the unwanted oscillations from the discrete evaluation, /, is smoothed by 
adding a second derivative term, i.e., /, = 75/,+. 125(/, +1 -/,_,). By default, two smoothing passes 
are made, but this can be overridden by changing the filter input control parameter, NFILT. The 
tension parameter, co,, is smoothed in the same way. 


1.3. 1.1 A and B. As seen in Eq. (3), co is also a function of A and B. These variables 
are the important elements of the self-adaptive nature of the scheme because they control the 
size limits of the grid spacing. A and B are computed from the user-supplied maximum and 
minimum allowable grid spacings, As max and A s MIN . These are relative values of mesh size and 
allow the user to specify the proportion of largest to smallest As in the final adapted grid. 

The value of A is constant throughout the grid adaption and is given by 

A = - 1 (10) 

MIN 

The value of B is computed (by an iterative process) for each j line to provide 

input As WN = computed As min 

That is, the computed minimum grid spacing is equal to the user-requested value. Appendix I 
gives a detailed description of the derivation of these two parameters. 

1.3. 1.2 Torsional effect on the evaluation of co. As just noted, the function of 
variables A and B is to control the size limits of the adapted-grid mesh spacing and, as described 
in appendix I, B is computed from the 1-D Eq. (2). Thus the addition of the torsion terms is 
somewhat inconsistent with the value of B, and the resulting grid spacings may give values 
slightly outside the requested limits. To provide for additional control, a weighting factor is 
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applied to co such that co, = (1 + A/®)co,., where co, is the weighted tightness factor. When equation 
(8) is solved, each As, is tested and if it does not lie within the requested spacing limits, co ([ will be 
computed; otherwise to, =1.0. If As, < As MW , we need to relax (decrease) the value of to, to 
produce a larger As, in the next iteration. If As, is too large, co, should be increased. We therefore 
use the modifier 


1. As, 

co, = -(- — — 
2 A s MIN 


1 . As, 


•+u 


MAX 


depending on whether As, < As MIN or As, > A s MAX . Equation (8) is therefore solved using the 
modified values of co. It should be noted that co, = 1.0 in the regions of side-spacing control (see 
section 1.4.3) where As is permitted to be outside the specified limits. 

1.3.2 Torsion terms. As shown in Eq. (7), the torsion terms, t,(s'-s,) and \|/,(s, -s { ) are 
added to the 1-D equation to preserve the required smoothness and orthogonality of the grid. 
They are constructed as being analogous to the behavior of torsion springs. Consider first the 
torsion term within the plane, f(x). A torsion force, T, relates the node at (i,j-l,k) to the node 
(i j k), is proportional to the turning angle, 0, shown in Fig. 3, and is defined as 

Ti * K0,.i-u < n > 


where k is a torsion-related constant. The next step is to relate this torsional force to the 
previously defined parameters on the adaption line. Since 0 is generally small, we can 
approximate 0, by (s' -s,)/\DA'\ (see Fig. 3), where (s'-s,) is a function of the inclination angle 
computed from the proportioned orthogonality and straightness vectors, as described in section 
I.3.2.3. In addition, we assume that K is proportional to the maximum co, (=1+A) and the local 
aspect ratio of the grid cell. Hence we can write 

_ fci + U-l (12) 

2|DA'| 


where X is the proportionality coefficient and is an input quantity. The variable t used in Eq. (8) 
can thus be defined as 


2\DA'f 


(13) 


The evaluation of \| f is similar, with the node (i,j-l,k) replaced by (i,j,k-l) and a 
proportionality coefficient X defined that is analogous to X. such that. 

_ ~ (14) 

V ' " 2\BA'f 

where |flA'| is the distance from (i,j,k-l) to the location of s’ jk . The evaluations of \DA'\ and \BA'\ 
are given in appendix II. 

1.3.2.1 Evaluation of s' and j* . These two variables, seen on the right-hand side of 
Eq. (8), are calculated from the torsion vectors and provide the direction of the smoothness and 
orthogonality constraints. The location of s' lies at the intersection of the projection of the 
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within-plane torsion vector t, with the arclength vector s, along the current adaption line. The 
location of s ■ is similarly determined from the between-plane torsion vector t*. Figure 4 shows 
the two unit torsion vectors: /, acting from point D and intersecting the line AC at A', and /*, 
acting from point B and intersecting the same line at A‘, such that 

s' = Sj + AA' and j* = .y,+AA* (15) 

In general, neither t nor F* will directly intersect the line AC, and the projection of the torsion 
vectors onto s is required. An example of the projection vectors is shown in Fig. 5 where the 
projection of F onto the k plane intersects s at A'. A full description of the calculation of the 
intersection is given in appendix II. 



Figure 4. Torsion vectors, t and V. Figure 5. Calculation of A', intersection 

of torsion and streamzvise vectors. 

1.3.2.2 Torsion vectors F and F*. Figure 6 shows in detail the within-plane torsion 
vector F, and the associated base vectors (n = f(u,b), and e). The vectors in these figures are 
represented on a 2-D surface, and it should be remembered that the shown f is actually the 

projection of f seen in Fig. 4. The unit vector t , associated with the nodal point (i,j,k) [but acting 
from the point (i,j-l,k)] is defined as 

+ (16) 

where 

C,is the input parameter that defines the percentage of straightness to orthogonality, 
e, is an average straightness vector from (i,j-2,k) to (i,j-l,k), and 

ii, represents the orthogonality vector between the j-1 line and the node (i,j,k) and is a 
function of b and u vectors described below. 
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The between-plane torsion vector 7* which 
acts from the node (i,j,k-l), is similarly defined as 


l i,j.k 




(17) 


where C* is the input proportion parameter, and 
£' and n* are the equivalent straightness and 
normal vectors. 

It should be noted that Figs. 4 and 6 show 
vectors intersecting the j line in the interval 
(i,i+l); however, this is not necessarily the case. In 
general, we assume that the projection of t, 
intersects the j line in the interval (1,1 + 1). 
Appendix II describes a general method used to 
compute the intersection of a vector with the 
vector s. 



1.3.2.3 Evaluation of vectors s, e, and n. The unit vector I, (direction AC in Fig. 7) is 
the arclength vector from s, to J i+1 along a j line. It is defined as 

kj= s J + s yj + s zk 

where 

s x = 5 y. an< l s z, = ^( Z i+lj.k~ Z i.j.k) 

The unit vector e, (direction ED in Fig. 7) is defined as the sum of three straightness 

vectors, d,_ t , d i and d M , where 

d t] =d x i +d y j + d z k 

such that 

d x< =n-\(Xi.j-ik- x i,j-2,k)' d y. =-^(y,.j-ik-yi.j-2.k) 


A(i,j,k)^ 



□ k plane 
k-2 plane 


and 


Therefore, 


d. ^i.j-2.k ) 


€i — ( d j-\ + d, + d j+ j ) 


If the line j-2 does not exist (such as at the 

Figure 7. Straightness vectors, e and e . second line from a physical boundary), it is 

assumed that e = h . 


8 






The straightness vector between planes, e *, is shown as PB in Fig. 7. It is computed like e 
above, with points (i,j-l,k) and (i,j-2,k) replaced with (i,j,k-l) and (i,j,k-2) respectively. 

The vector h is a combination of two vectors: 



"=nM+a-'j«] as) 
W l J 

where u, represents orthogonality to the ik 

plane through the j-1 line, and b t the 
orthogonality to the ik plane through the j 
line, as shown in Fig. 8. The parameter t n 

proportions b and u and is defaulted to 0.5, 
changing only at the upper and lower 
marching boundaries, as discussed in section 
1.4.4. The calculation of a normal vector to a 
point in a 3-D surface needs defining: for 
example, in Fig. 8, the vector u acts at D, 
perpendicular to the ik surface at j-1 (the 
surface that contains G, F, M, and Q). In this 
case, there are four possible planes to choose 
from (MDF, MDQ, GDF, and GDQ) and thus the 


definition of the u and b vectors is the 
average of the normals to each of the four planes, when the planes exist. Since two vectors define 
the plane in which they lie, the normal to that plane is obtained by taking their cross product. 


The vector b would also act at D but is perpendicular to the ik plane at j. Appendix III describes 
the general calculation of the normal to a given plane at a given point. 

In addition to the within-plane orthogonal constraints, orthogonality is also controlled 

between planes by a function of the vector n , defined as 




(19) 


The normal u acts at point B, perpendicular to the ij plane at k-1, and b acts at A, but 
perpendicular to the plane at k. 


1.4 Enforcement of Boundary Conditions. 

The ability to modify the adaption techniques in boundary regions substantially improves 
the flexibility of the adaptive scheme. The adaption domain is defined by the user and may be 
equal to, or a subset of, the physical grid. A boundary occurs at the limits of the adaption domain 
defined by the user. Within a plane, there are two types of boundaries: marching boundaries (all 
points along the initial and final adaption lines) and edge boundaries (at i = i st and i = i end )- There 
are also initial and final surface boundaries to take into consideration. Figure 9 shows a 2-D 
example of an adaption domain as a subset of the physical grid, and illustrates the two types of 
boundary lines when marching in the j direction. Note that if marching had been in the i 
direction, the boundary definitions would have been reversed. 
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By internally amending the previously defined 
variables of C,, t n , and X in the boundary regions, it is 
possible to maintain physical characteristics at 
boundaries, allow for multiple passes, provide 
continuity across grid boundaries, and simplify 
multiple grid adaptions. Specialized boundary 
conditions, such as wall shape preservation and 
periodic boundaries, are discussed in separate 
sections. 

1.4.1 Treatment of initial marching boundary 
line. If the initial adaption line is within the physical 
domain, a smooth transition will be required across 
the starting internal boundary. For example, this 
situation occurs when adapting in zones; each zone 
has different flow features and the user may wish to 
march up to a certain line using one set of 
parameters, then continue marching using a new set. 
The common boundary between the two zones must 
remain unchanged when starting the second adaption pass. This feature is controlled by the 
input parameter m g . When m g >0, a smooth transition from external grid lines (e.g., those in the 
already adapted zone) to internal lines is created by maintaining the same grid points along the 
initial line and then incrementally introducing the input adaption parameters. For example, 
when stepping to the second adaption line, as much straightness as possible is maintained by 
setting C, = 1 (see Eq. (16)). At each subsequent line, C, is gradually decreased until it equals the 
user-supplied input value. The number of lines stepped before the full effect of the new 
parameters is felt depends on m g ; C, will linearly decrease until nt g lines have been adapted. An 
example of this can be seen in Fig. 9, in which j sl =4 and =5. The grid points along j-A will 
remain unchanged while those along j — 8 will be fully adapted to the input control parameters. 
The code controls the adaption parameters for lines in between. Consequently, we define a 
variable n m as 



m -j 


and replace the value of C, used in Eq. (16) by 


C, = C,(\-n m ) + n, r 


( 20 ) 


At the same time, the value of X is increased to >.(l + 5nj (thus increasing the magnitude of the 
torsion term) so that this amendment to C, is effective. After m g lines, n m = 0, thus returning 
C,and X to their original values. Note that the computation of e along the j s , + 1 line is a function 
of the j st - 1 line (if it exists), and this will also help in the merging process. 

A similar parameter, m' g , is used when the first plane must remain unchanged (e.g., at a 
multigrid boundary) and subsequent planes are gradually adapted. 

1.4.2 Treatment of final boundary line. Another discontinuity will occur between the final 
adaption line and any subsequent lines external to the adaption domain. Since the adaption 
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process is a marching scheme, it is not possible to use the same merging concept described earlier. 
On request (through the MARCH input parameter), the grid points on the remaining external 
lines will be redistributed with the same proportions as the points on the final adapted line. In 
addition, nonadapted planes can be proportioned with respect to lines on previously adapted 
planes. This is not an adaption to the flow field, but provides a more acceptable interface between 
the computational and physical domains. 

1.4.3 Treatment of side-edge boundary. A smooth transition will also be needed at the 
side-edge boundaries if the adaption domain is internal to the given grid, i.e., if i st > 1 or 
i end <imax. Figure 10(a) shows an example in which the fixed external grid spacing is denser than 
the first redistributed points, giving a discontinuous effect across the boundary. If the user 
requests continuity of mesh spacing across the side-edge boundaries (by setting the input 
parameter NEDGE, (n g )*0), a modification is made to the tension parameter, a>. Consider the 
start-edge boundary at i = i sl along line j. We wish to enforce some value on co, that will give a 
value of Aj, close to the value of Aj ( ,. To do this, we find the average coAj along the converged 

j-1 line and replace to, by coAj/Aj, _ I; . 

This value remains fixed during the 
iteration, but is merged into the updated 
values of to close to the boundary. In 
general, this implies that co 2 , co 3 and co 4 
will be amended, however an additional 
option is available (using input 
parameters MG1 and MG2) to spread the 
effect further into the adaption domain. 

Figure 10(b) shows the effect of this 
process, giving a more appropriate 
spacing in the vicinity of the i sl 
boundary. The end-edge boundary is 
handled in a similar manner. 

The side-edge spacing often needs 
to be improved even when the adaption domain coincides with the physical domain. If an 
adaption pass generates inappropriate side-edge spacing, the code can be rerun with n g set to non 
zero in an attempt to improve the spacing by using the above technique. In this case we do not 
have an external As, and Aj 2; is used instead of Ay, • This will usually prevent the spring 
constants from pulling the lines too far off the boundary. 

The variable n g can be set to initiate computation for either or both side edges. 

1.4.4 Treatment of orthogonality at boundaries. The code provides the choice of 
constructing grid lines that are as orthogonal as possible to a marching boundary, either to the 
final line in a plane or to the entire final plane. To accomplish this, the normal vector ft (and/or 
n *) is emphasized over the straightness vector (by decreasing C, when the adaption line is close 
to a marching boundary line. For even greater control, the coefficient t n is modified to emphasize 

either u or b, depending on whether the initial-line or end-line boundary is being considered. To 
ensure that the modifications to these torsion coefficients sufficiently affect the computation, the 
value of X is simultaneously increased to accentuate the torsion term. Since not all marching 


side edge boundary 


side edge boundary 


O 1 n~\ 

(a) ^ (b) S 
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Figure 10. Control of side-edge boundaries. 

(a) with no edge control; (b) with edge control. 
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boundaries are physical boundaries, an input option is available that will override this emphasis 
on orthogonality for either boundary. 

1.5 Preservation of Initial Wall Boundary Shape 


In applications for which the shape of the wall boundary is defined by large geometric 
gradients, such as for turbine blades, sharp corners, and the leading edges of airfoils and wings, 
sufficient points need to be placed in the appropriate regions of the initial grid to accurately 
define the geometry. If flow-field gradients are weak in these regions, the standard grid- 
redistribution algorithm will cause points to be dispersed, leaving insufficient points to properly 
describe the boundary. To maintain necessary clustering in these regions, on user request a new 
variable is introduced that is a function of the geometry gradients that define the boundary shape. 
The weighting parameter now becomes a function of both flow-field gradients and geometry 
gradients. The solution procedure will therefore redistribute the points into regions of hig 
geometry gradients as well as into regions of high flow-field gradients. Both the start and end 
wall boundaries are treated. A merging technique is used to integrate the boundary and internal 
redistribution, so that the internal flow is controlled only by the flow-field gradients. 

The original weighting function was defined in Eq. (3) as 

co= 1 + Af B 

where f is a function of dq/ds. When the geometry option is invoked, / becomes a function of 
both dq/ds and dg/ds, where f(dg/ds) is computed as the radius of curvature at the wall 
boundaries. The radius of curvature is defined as 


dx d 2 y dy D 2 * ^ j ( 3* 2 + dy 1 ^ 

35 ds 2 d 5 35 2 J/ ds d s ) 


where the derivatives are computed from the spline coefficients. 

The full effect of the geometry function should be felt at the wall boundaries; it should not 
be a factor in the internal grid redistribution. Hence the final weighting function takes the form 

co = 1 + 

where f = C q f(q) + C g f(g) and the constant C q normally equals 1. The constant C g equals 1 when 
the upper and lower wall boundaries are being adapted, but it must be decreased away from the 
walls until C g = 0 internally. The value of the aspect ratio was chosen to drive the rate of decrease 
of C; i.e., C = \-AA s , where A s is the aspect ratio at the maximum radius of curvature. 
Regardless of the value of A s , C g remains positive for a minimum of four steps from the 
boundaries and decreases smoothly. When the upper boundary is approached, C g must be turned 
on when necessary and increased with each step, and f(g) becomes a function of the upper wall 
boundary only. 

It is also possible to adapt only to the geometry gradient. If this is requested, C q — 0 and 
C =1 throughout the flow, and f(g) is computed for each adaption line, based on the local 
geometry of that line. This option has proven to be a useful tool for improving the starting grid 
before any solution is obtained. Points will be smoothly clustered with respect to the geometry 
gradients. 
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1.6 Finite-Volume Grids 


In applications using finite-volume techniques, each flow-field variable (q) is evaluated at 
the cell center instead of the nodal point of the grid. These cell centers are usually positioned at 
the average of the surrounding grid points (four for 2-D and eight for 3-D). This creates a solution 
file that has one less point in each direction than the grid file. The finite volume module in 
SAGE interpolates the q values onto the grid points before performing the adaption. On 
conclusion the q values are re-interpolated at the new cell centers. 

A feature of finite volume grids is the definition of the outermost boundary cells in all 
computational directions as 'ghost cells' containing specified flow values. The adaption 
procedure destroys these values and SAGE employs a very simple method to replace them: on 
output, each ghost cell contains the same value of q as the adjacent internal cell. It is therefore 
highly recommended that the user re-apply the boundary conditions before processing the new 
grid. This is true for both 2-D and 3-D grids and implies that if a user wishes to adapt one plane 
for testing purposes, an internal plane would be more appropriate. 

1.7 Multiple Grids 

When investigating the flow around complex structures, computational grids are fre- 
quently organized in multiple-grid format. This enables each individual grid to be of manageable 
size while maintaining a single dataset and also allows for overlapping grids. The original SAGE 
code could only adapt single grids and it was therefore necessary to separate the grids before 
adaption. Version 2 of SAGE can read and write datasets in multigrid format and also provides a 
new input control parameter to specify which grid to adapt. In addition, a major new feature in 
the SAGE code is the capability to transfer data between matching zonal boundaries (with 1:1 
mapping) in separate grids. Some or all boundaries of each individual grid will match in some 
part to boundaries in one or more of the other grids. It is important that common boundaries, 
where data in separate grids represent the same location in computational space, retain the same 
grid and flow distributions. This can be seen in the simple two-part multiple grid shown in Fig. 
11(a), where the shaded plane is common to both grids. After adaption, both planes must still 
contain matching data. 



Figure 11. Simple multi-zoned grid, (a) Common plane; (b) restricted adaption direction with 
original SAGE; (c) preferred adaption direction, using plane transfer procedure. 


The adaption technique used in SAGE is a marching scheme, and therefore the order and 
direction of adaption have a marked effect on the final grid-point redistribution. The only way 
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the original SAGE code could adapt the two grids and still retain the common boundary data is 
shown in Fig. 11(b): the common plane had to be adapted first, with plane marching occurring in 
opposite directions. The new export/import feature in version 2 enables the order of adaption 
shown in Fig. 11(c): after grid 1 has been adapted, adapted data from the final plane can be 
transferred to the first plane in grid 2. Then, using the merging feature (MGPLS) to prevent the 
first plane of grid 2 from being corrupted, adaption can continue in the same plane marching 
direction. All three steps (adaption of grid 1, data transfer, adaption of grid 2) can be accomplished 
in the same execution pass. Figure ll(c )shows only the simplest multigrid structure. In reality, 
multiple grids may contain many grids, with partial matching surfaces to more than one grid. 
Input options are available to handle this and examples can be found in Section 3. 

1.8 Discussion 

1.8.1 Flexibility and segmentation. The vector approach used in the analysis provides for 
great flexibility in the use of the SAGE code. The user can choose the adaption direction and order 
of sequential adaptions without concern for the computational data structure. Multiple passes are 
available with no restraint on stepping directions: for each adaptive pass (or sweep), the user can 
choose a completely new set of adaptive parameters. This facility, combined with the capability of 
edge-boundary control, enables the code to handle multidimensional zonal grids without losing 
continuity along the common boundaries. For patched grids, the multiple-pass capability allows 
complete adaption. It was mentioned in the introduction that the adaptive technique does not 
produce a unique grid. This non uniqueness enables the user to choose the most appropriate 
solution to the grid adaption. 

1.8.2 Blanked grids. Certain types of complex multiple-grid structures utilize the facility 
in PLOT3D that blanks out regions of overlapping grids. When this option is invoked, the 
number of points in each grid line need not be constant, i.e., IMAX can vary for each line. A 2-D 
version of SAGE has been modified to handle this situation, but in most applications it is not 
practical. The major drawback occurs when the number of points on a line decreases with each 
step. Depending on the flow-field gradients and the currently adapted line, the "dropped" points 
can produce a discontinuous grid structure. The starting location of the adaption can be chosen to 
ensure only an increasing point count, but this removes flexibility. This blanking feature is 
therefore not included in the standard version of SAGE since it detracts from the simplicity of the 
code. Nevertheless, the blanking option has become a useful tool for complex overlaid grid 
structures, and it is anticipated that the next version of SAGE will contain this option. It is 
suggested that the authors be contacted for further information. 

1.8.3 Cyclical and Periodic Boundaries. Certain classes of grids used by flow-field solvers 
are not suitable for adaption using the basic grid formulation described thus far. The self-adaptive 
grid procedure is a marching scheme, i.e., the solution along each line is influenced only by 
previously adapted lines. The adaption of the first grid line is based on flow gradients only, but as 
marching proceeds, the redistribution of points is dampened by the torsion effect. As a result, the 
final adaption line will be somewhat less adapted to the flow-field gradient than the initial line. 
This implies that the scheme cannot be directly applied to cyclical and periodic grid structures, 
since common and/or matching boundaries will show as discontinuities. For 2-D problems, these 
difficulties have been overcome by including both an iterative and a rearrangement procedure 
that can be requested by the user. However, these options are not provided in this version of 
SAGE even when TWOD is specified, and interested users should contact the authors to obtain a 
copy of SAGE2D and its documentation. For 3-D grid structures, maintaining periodicity at 
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matching planes is a more difficult task. Since the initial grid and flow-field on matching first 
and last planes are identical, setting the plane torsion parameter (X ) to zero produces matching 
adapted grids on these planes. In this case, the adaption of intermediary planes has no influence 
on subsequent planes. 


1.9 Appendices 


1.9.1 Appendix I: Derivation of A and B. 

A and B provide self-adaptiveness to the adaption scheme. These two parameters ensure 
that all repositioned grid nodes maintain grid spacing to within the user-requested mesh-size 
limits (A s M/N and As max ). 

1.9.1.1 Calculation of A. We wish to relate A to the input values of A s MIN and 
A s MAX . A is constant throughout the entire mesh, and hence the 1-D relationship given in Eq. (2) 
holds. From the original definition of / in Eq. (9), f max = 1 and f mm =0; therefore, from Eq. (3) we 
have © ma * = l + A and © mm = 1 . From Eq. (2) (rewritten as A j, = KYco ( ) we can see that the 
minimum As will occur at K/at^. Similarly, the maximum As occurs at K/(o mm . We therefore 
wish to set A s MIN = K/(\ + A) and As max = K. These can be solved simultaneously: by eliminating 
K we get the expression given in Eq. (10) 

A _ ^max t 
MIN 

1.9.1.2 Calculation of B. This calculation is more complex. B is found by an iterative 
procedure and will change for each j line. The objective is to determine which value of B will 
give the minimum computed As, equal to the requested minimum, A s MIN . For each value of B 
there exists a minimum As,. (An example of a plot of B vs. As mi „ is shown in Fig. 12). We need to 
find B at the requested A s MlN . To do this, we assume an initial value of B (= 1.0), evaluate ©, and 
then solve for the new As using Eq. (5). Although Eq. (5) is true only for the initial line, it is 
assumed here to hold for all j in order to simplify the B calculation. If |As MW -minAs{ ,, '| is small, 

an acceptable value of B has been found. If not, a new B is computed from B* n+l> = B (n> + AB (n> , co is 
reevaluated, and the procedure repeated. 

B <n> can be found from the definition 


— min( As,) 


lim (& s min ~ minAs[ n) ) 
dB-> o . A B (n> 


( 21 ) 


As mentioned in the calculation of A, we know that As, is a minimum when co, =1 + A, 
and by substituting this in Eq. (5) and differentiating, we obtain 


where 


3 e 

— min(As.) = — WL 

dB ' 1 + AdB 


aro 


<p 


( 22 ) 
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/A 
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We know that 


_d_f J_1 1 3<p 

3BV<p J <p 2 3 B 

Hence, the next step is to evaluate 3<p/3 B. With this definition of (p, we can take the summation 



from Eq. (5), we obtain 


—min(As,) = 


A(\ + A) 




|2 K'fljogf, 


i=i 


(a. 


Finally A B can be obtained from Eq. (21), giving B. 

1.9.2 Appendix II: The intersection of a 3-D vector, v, with the arclength vector, s . This 
technique is used to obtain the b vectors used in the orthogonality terms and to find the location 
of s' and s*. For example, to determine s' we need the segment in which the torsion vector t 

from the grid point (i,j-l,k) intersects the j line. To evaluate the direction cosines of the b vector, 
we need to find the segment / — > / + 1 along the j line that contains the intersection of a vector 
acting from (i,j-l,k) that is normal to the; segment. For s' and b* , we are concerned with the line 
) in the plane k-1. In all cases, the technique for finding / is the same. 

Figure 5 shows the case of the torsion vector f, computed for the point (i,j,k) and 
emanating from the point (i,j-l,k) at D. This vector (/) most likely will not lie in the plane 
described by ADC (or the equivalent 1 plane) and thus its projection onto the plane is required. 
The value of s' in the solution equation lies at A ' , the intersection of the projection of t onto the 
ADC plane and the arclength vector, s. The vector n (=n x i +rt y j + n z k) is the unit normal to the 
plane ADC and is computed by the technique described in Appendix III. Both |A4'|and \DA'\ are 
required in the code. 

From vector addition, we can write 


|AA'|f = AD+ \DN\t - \NA'\b 


(23) 
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Since the coordinates of A and D are known, AD(=a x i + a y j + a z k) can be evaluated, hence 

\AA'\ |0/V| and |M4'[ are the three unknowns. By equating coefficients of i, j, and k in the vector 
Eq. (23), we can obtain a set of three equations with three unknowns: 

| AA'\s x =a x + |DAf|r x - |M4'|n x , | AA'\s y = a y + \DN\t y - \NA'\n y and \AA'\s z = a z + |ZW|f, - 1 NA'\n z 

These equations can be readily solved for \AA'\ by computing the relevant determinants. 

\DA’\ is also found by utilizing vector addition. We have 

DA' (= \DA'\DA ') = \AA'\s - AD 
and, since the lengths of the left- and right-hand sides are equal, 

\DAf = (\AA'\s x -a x ) 2 + (\AA’\s y -a y ) 2 + (\AA'\s z - a z ) 2 

1.9.3 Appendix III: Definition of normal vectors. The vector normal to a plane is required 
for the calculation of vectors u, u', b, and b ’ . This appendix describes the derivation of a general 
unit-vector normal, n, at a given point (i,j,k) and normal to the plane ik passing through the 
constant j-1 line. This specific orientation is analogous to u. Figure 8 shows such a vector acting 
at D, that represents the normal to the shaded surface. In the figure it can be seen that there are 
four possible planes passing through (i,j-l,k): GDF, FDM, MDQ, and QDG. Any vector normal 
required in the code is calculated as the average of the normals to each of these four planes. 
(When the node (i,j-l,k) is on or close to a boundary, only one or two of these planes will exist.) 

The normal to a plane is derived from the cross product of two vectors lying within the 
plane, with care taken to ensure the correct order of the operation (this is why the input grid 
should be organized as a right-handed system). In the example shown in Fig. 8, the four cross 

products are DFxDM, DFxGD, QDxGD, and QDxDM. The required unit normal is thus the 
unit vector representing the sum of these four vectors. 

This picture gives an idealized view of the defined planes. In reality, some of the required 
lines have already been adapted and some have not. In Fig. 8, the point F could be significantly 
different from D, since D is an already-adapted point and F is still part of the initial grid. The code 
solves this problem by forming a block of data around the current j line that includes all i for 
lines j — 1 => j + 1 for planes k — 1 => k + 1. Lines that have not yet been adapted (e.g., j+1 at k or j-1 
at k+1) are proportioned with respect to the s, at j; this relocates the points for a smoother 
computational effect. 
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2. SAGE USER GUIDE 


2.1 Overview 

The SAGE code is based on the self-adaptive grid method developed by Nakahashi and 
Deiwert (1985). This guide contains information that will enable the user to run the code with 
little knowledge of the mathematical concepts employed in the development of the adaption 
technique. Included is a detailed description of the input-control parameters, along with 
routine descriptions and nomenclature. The code stands alone and has been run on many 
computer systems. Users wishing to understand and/or amend the code can find the details of 
the mathematical background in the first section of this document. 

The first step in the code is to read two data files: one that contains the coordinates of 
the grid (x,y,z) and another that contains the corresponding flow-field variables, cj. It is 
assumed that both these files are in the binary format (given on the next page) associated with 
the plotting software package PLOT3D. A single-grid, finite-difference, 3-D format is assumed 
unless re-specified by the user. The next step is to adapt the grid with respect to the user- 
supplied input-control parameters. Adaption takes place as a sequence of one-directional 
adaptions, with the input-control parameters determining the order of adaption. Each plane 
(i.e., 3-D surface) is adapted in one direction only before stepping to the next plane. Figure 1 
(see* page 2) shows a small section of a 3-D grid: adaption in the i direction (stepping in j) has 
already been performed on the first plane. The code then steps to the next plane and performs 
the same directional adaption on this plane. A 3-D adaption is performed by a sequence of 
adaptions: once the grid has been adapted in one direction, the adaption can be repeated with a 
new parameter-input set describing another direction. (It should be noted that different orders 
of the adaption direction will produce different adapted-grid solutions). When all adaptions 
are complete, the code sends the redistributed grid points and the corresponding interpolated 
flow variables to two new data files, also in PLOT3D format. 

The analysis that generates the algorithms used in the code is given in detail in the first 
section of this report. Briefly, the redistribution of points is controlled by parameters related to 
torsion and tension springs, and by the maximum and minimum allowable grid spacings. 
Along each grid line, (x,y,z) is transformed to a 1-D arc-length variable, s. If we assume that 
adaption (within a plane, k) steps in the j direction, the tension spring constants, co, are 
evaluated at each point (i,j,k); i.e., co, = f(s, j k ) and are a function of the gradients of the user- 
chosen flow-field variables. The torsion parameters (shown in Fig. 2) x ( - = f(s l j k ,s i j _ lk ,s i j_ 2 ^) 
and \jc, =f(s, jk ,s l)k _ l ,s l j k _ 2 ) are the link between the current line and the previously adapted 
lines. These are the parameters that maintain the integrity of the grid by controlling 
straightness and orthogonality within the plane and between planes. With the evaluation of 
these variables, a system of («, -2) equations (given in Eq. 8) is developed for the current j line 
and solved as a tridiagonal system. Once the adaption is complete for the current j line, the 
code steps to the next j line (either forward or backward) and repeats the same procedure. At 
the end of the current k plane, the code moves to the next plane and the process is repeated. 

There are eight major steps taken in the code: 

1. Input of three data files: initial grid, flow-field solution, and user-control parameters 

2. Initialization and reorganization of data for computational purposes 

3. Adaption along the start line on the start plane with the 1-D technique 
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4. For each subsequent j line on the initial plane; 

a. Computation of variables that define the torsion and tension springs (co and x), 
and hence coefficients of the tridiagonal matrix defined by Eq. (8) 

b. iteration to find new values of s, and hence (x,y,z) 

5. Repetition of step 4 until all lines are complete on this plane 

6. For each subsequent k plane, computation of the torsion spring \|/ for inclusion in 
the coefficients of the solution Eq. (8) 

7. Repetition of steps 4 through 6 until all planes have been adapted 

8. Output of new grid and interpolated flow-field files in the original format 

2.2 Execution of the SAGE Code 


The following is an example of the command file to run the SAGE code on a UNIX 
system. SAGE is written in FORTRAN and is self-contained; sage is the executable module. 

cp xyz.gr d fort. 7 ! copy grid file to unit 7 

! copy solution file to unit 8 

! run SAGE code with sage.inp containing user control parameters 
! name the output adapted grid file 
! name the interpolated function file 


cp q.futi fort. 8 
sage < sage.inp 
cp fort. 10 xyz.out 
cp fort. 11 q.out 


where sage.inp is in namelist format ($NAMEL). 


{Note: Users on Silicon Graphics Inc. machines will need to amend the READ(5,NAMEL,...) 
statement in subroutine INITIAL to READ(5,NML=NAMEL / ...)J. 


The remaining files are in PLOT3D format: 

xyz.grd contains the initial grid points (stored as a right-handed coordinate system) 
q.fun contains the flow-field variables to which the grid is to be adapted 
xyz.out contains the adapted grid points 

q.out contains the flow-field variables interpolated on the adapted grid 

In addition, an output message file is associated with unit 6, which is often defaulted to the 
user's output device. For other operating systems, the user must appropriately assign the six 
input/output files. 

The following are the read statements for single-grid 3-D PLOT3D binary/unformatted 


input files: 

xyz.grd: 


READ(7) 

IMAX,JMAX,KMAX 

READ(7) 

(((X(I,J,K),I= 1 ,IM AX) ,J = 1,JMAX),K= 1,KMAX), 
((( Y (I,J,K),I= 1 ,IMAX),J = 1,JMAX),K= 1,KMAX), 
(((Z(I,J,K),I=1,IMAX),J=1,JMAX),K=1,KMAX) 

q.fun: 


READ(8) 

IMAX,JMAX,KMAX 

READ(8) 

FSMACH,ALF,RE,TIME 

READ(8) 

((((Q(I,J,K,N),I=1,IMAX),J=1,JMAX),K=1,KMAX),N=1,NDIM) 

where 



IMAX=number of points in the i direction of the grid file 
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JMAX=number of points in the j direction 
KMAX=number of points in the k direction 
NDIM=number of Q variables (default=5) 

As seen in the above format, five flow-field variables are expected in the Q file. PLOT3D 
preassigns p, p u, pv, pw and e , but since the SAGE code requests only the index of the 
function, any variables may be stored. Note that it is possible to handle any number of flow- 
field variables by changing the value of NDIM in the parameter statement at the beginning of 
each subroutine and recompiling. Also contained in the parameter statement are the grid 
dimensional variables ID, JD, KD and IMX, currently set at 75. These are used to define the 
required maximum dimension of the grid arrays. Because of the internal switching of data, 
these values may not exactly coincide with IMAX, JMAX, and KMAX. If the assigned code 
dimensions are too small, a message will be sent to the user stating the minimum dimension 
requirements for the current application. If the grid is in multiple grid format, an additional 
record occurs at the beginning of each file stating the number of grids and the grid dimensions 
are supplied in dimensioned arrays. FSMACH, ALP, RE, and TIME are not used in the code 
and may contain dummy values. They are part of the PLOT3D package and are displayed on 
the output plots. 


2.3 User Input Parameters 

The file sage.inp is the user-supplied, input-parameter control file. The grid adaption is 
based on the user's choice of these input parameters, which are listed and briefly described 
below. This is followed by a more complete explanation of each parameter. The input file 
sage.inp uses the namelist format, since generally only a few of the input parameters need to 
be changed from the default value set by the code. These default values are shown in 
parentheses in the list given in section 2.3.1 and must be input by those using a nonstandard 
version of Fortran. If more than one adaption pass is to be made (for example, a two- or three- 
directional adaption), the subsequent adaptions can be made on the adapted grid by linking up 
to 10 sets of namelist inputs within the same sage.inp file. For multiple grids, these sets can 
also contain multiple export/import passes. 

Before describing the input parameters, some terminology needs to be clarified: 

The term physical domain is used to reference the complete grid defined by the input 
grid file, i.e., the grid bounded by IMAX, JMAX, and KMAX. The adaption domain is the part 
of the grid, as defined by the input-control file, that is to be adapted, i.e., (IST,IEND), 
(JSTJEND), (KST,KEND). These two domains can be equivalent. The direction h j, or k refers 
to the direction of the grid coordinates as defined by the order in which they are stored in the 
grid file. The first index (normally containing x) is named the i direction, the second index is 
the j direction, and the third index the k direction. This implies that if data happened to be 
stored as (z ,x,y) instead of ( x,y,z ), i would represent the z direction. Regardless of the order, the 
data must be stored as a ri ght-handed coordinate system. The adaption direction is used to 
define the 1-D line along which adaption (redistribution of points) takes place. In the analysis 
and in the descriptions below, this direction is always i for convenience, but the user may 
request i, j, or k. There are two stepping directions : one within the plane, defined as j in the 
analysis, and the other in the direction of plane marching, defined as k. Although the analysis 
and descriptions in this report assume this particular order of adaption and stepping, the code 
makes no such a priori assumptions, since the order is controlled by the user's input- 
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parameter file. Reference is also made to grid boundaries . Side-edge boundaries refer to the 
start and end points at the edge of the adaption line (assumed to be i in the analysis). These are 
the edges that are controlled by the NEDGE parameter. Marching boundaries are the start and 
end lines of the stepping (within the plane) direction (assumed to be j). Adaption close to these 
lines are affected by the parameters ORTHS(l), ORTHE(l), and MGSTEPS, as well as by some 
internal controls. Planes juxtaposed to start and end planes are affected by ORTHS(2), 
ORTHE(2), and MGPLS. 

2.3.1 Parameter control file, sage.inp 

The input parameters, with their default values (in parentheses) and short descriptions, 


are 



$NAMEL 

1ST 

(1) 

first adaption line in i direction 

IEND 

(IMAX) 

last adaption line in i direction 

JST 

(1) 

first adaption line in j direction 

JEND 

(JMAX) 

last adaption line in j direction 

KST 

(1) 

first adaption line in k direction 

KEND 

(KMAX) 

last adaption line in k direction 

IJPLANE 

(TRUE) 

the adaption surface lies in the (i,j) plane, with 
plane stepping occurring in the k direction 

JKPLANE 

(FALSE) 

the adaption surface lies in the (j,k) plane, with plane stepping 
occurring in the i direction 

IKPLANE 

(FALSE) 

the adaption surface lies in the (i,k) plane, with 
plane stepping occurring in the j direction 

ISTEP 

(FALSE) 

=.true. for stepping in the i direction within the plane: it cannot 
=.true. if JKPLANE =.true. 

JSTEP 

(TRUE) 

=.true. for stepping in the j direction within the plane 

KSTEP 

(FALSE) 

=.true. for stepping in the k direction within the plane 

INDQ 

(1) 

index of adaption flow- field variable, q; default of l=>p 

IQ(8) 

(0) 

enables a combination of q variables to drive the adaption 

RDSMAX 

(2.0) 

relative maximum allowed Aj:Ay M4X >1.0 

RDSMIN 

(-5) 

relative minimum allowed As:As MIN < 1.0 

CLAM(l) 

(.01) 

=X, coefficient of torsion parameter x, 1 < X < 10 -6 

CLAM(2) 

(.0001) 

= A.*, coefficient of plane torsion parameter, \|/; same range as X 

CT(1) 

(.5) 

proportion of "straightness" to "orthogonal" for torsion vector F 

CT(2) 

(.5) 

as CT(1), but proportions plane torsion vector, F‘ 

NEDGE 

(0) 

override control on side-edge adaption: 
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MG1 (0 or 4) used with NEDGE: number of points to merge start-side spacing 

MG2 (0 or 4) used with NEDGE: number of points to merge end-side spacing 

INTER (2) order of interpolation: 2, 3, or 4 

NFILT (2) number of passes to filter q and w 

MGSTEPS (0) number of merging lines for within-plane torsion controls 

MGPLS (0) number of merging planes for between-plane torsion controls 

MARCH (FALSE) =.true. to extrapolate end adaption line throughout remaining lines 

MARCHPL (FALSE) =.true. to extrapolate last adapted plane throughout remaining planes 

ADD (0) =n to add n points between each node in selected range 

LSTADD (1ST) lower limit of range for adding points (i.e., when ADD>0) 

LENDADD (IEND) upper limit (if ADD=0, values are ignored) 

SUB (0) =n to delete n points between each node in selected range 

LSTSUB (1ST) lower limit of range for point deletion (i.e., SUB>0) 

LENDSUB (IEND) upper limit (if SUB=0, values are ignored) 

REMOVE (0) removes requested number of points from outer grid region 

NOUP (FALSE) =.true. if no adaption required (e.g., for adding points only) 

SAVE (TRUE) =.false. to suppress output of data files 

ORTHS(2) (TRUE) =. false, to remove orthogonal constraint at start boundaries 

ORTHE(2) (TRUE) =. false, to remove orthogonal constraint at end boundaries 

GEOM (FALSE) =.true. to include geometry at wall boundaries in adaption variable 

QFUN (TRUE) =.false. to remove f(q) from adaption variable 

(used with GEOM=.true. only) 

NOQ (FALSE) =.true. if no q file is available 

LNSING (0) -n if this line is common to all adaption planes 

PLSING (0) =n if adaption plane n is collapsed to a line 

TWOD (FALSE) =.true. if datasets are 2-D (see section 2.4.1) 

FV (FALSE) =.true. if q file in finite-volume format 

MGRID (0) adaption grid number if multiple grids 

EXPORT (FALSE) denotes an export plane transfer, not an adaption pass, multi-grid only 

IMPORT (FALSE) denotes an import plane transfer, multi-grid only 

MPLANE (1) defines plane number, used for multiple grid export and import only 

IS,IE (0,0) defines i direction domain for plane transfers 

JSJE (0,0) defines j direction domain for plane transfers 

KS,KE (0,0) defines k direction domain for plane transfers 

First adaption Although many parameters have just been described, generally only a 
few are used for each adaption. The best technique, even for experienced users, is to retain as 
many of the default parameters as possible, view the results and then adjust some parameters 
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if necessary. For some 3-D problems, it is clear which plane should be the stepping plane so the 
PLANE parameter, along with the STEP parameter, can be initially chosen. Note that only one 
PLANE parameter and one STEP parameter are needed for one adaption pass. Also, a flow- 
field variable should be chosen for INDQ that shows the flow features most clearly. To find 
appropriate within-plane parameters, adapt on only one plane and/or turn off the 
connectivity between planes (i.e., CLAM(2)=0) to see the set of 2-D adaption planes. If there are 
any lines common to all planes or if a plane is collapsed to a line, PLSING and/or LNSING 
must be set. If the first adaption pass is not acceptable, try increasing the ratio between 
RDSMAX and RDSMIN, decreasing CLAM(l), and/or setting NEDGE=1. Since many of the 
parameters have interdependent effects, it is better to change only one or two of the 
parameters at a time. 

2.3.2 Explanation of user-supplied input parameters 

The following is a detailed explanation of each of the input parameters; they are listed 
in alphabetical order. 

ADD When this is set, points are added between adjacent mesh points within the 
requested range (see LSTADD, LENDADD). For example, if ADD=2, two points will be added 
between each consecutive grid point. Adding occurs only in the adaption direction and not in 
either of the stepping directions. Ensure that the added points do not cause the coded array 
dimensions to be exceeded. ADD and SUB can be used in the same pass to move points: note 
that the ADD will occur before the SUB, so the range on the SUB parameter needs care. 

CLAM CLAM(1)=A and CLAM(2)= A’ define the magnitude of the torsion parameters x 
and \|/, respectively. As the values of these parameters decrease, more points will be pulled into 
the high-gradient regions, at the possible expense of grid smoothness. 

CLAM(l) controls x, the torsion parameter within the plane; its order of magnitude can lie 
between KT 6 and 1. A value of zero produces a set of independently adapted lines, possibly 
generating crossed grid lines. As A. increases, the grid becomes smoother but less adapted. 

CLAM(2) controls the magnitude of vp, the torsion between planes. It has the same range of 
values as CLAM(l); however, if it is zero, the adapted grid may still be acceptable. For periodic 
planes (i.e., if the first and last planes are the same), setting CLAM(2)=0 is necessary to prevent 
a discontinuous grid at the juncture. 

CT CT(1)=C, and CT(2)= C*; they represent the direction of the torsion vectors Fand F* 
(whereas A and A* are their magnitude in these directions). They have the range of 
0< C,,C,* < 1.0, where a value of zero emphasizes orthogonality and a value of one emphasizes 
straightness. The default of .5 places the torsion vectors halfway between. This value is suitable 
in most cases, but it may cause problems when side boundaries are concave or when adapting 
on already adapted planes. 

EXPORT For multiple grids only. When EXPORT is set to .true., the user-supplied 
parameters for this adaption pass describe the current location of a plane of data to be 
transferred from one location to another, most probably in another grid. This parameter set 
must be followed by an IMPORT parameter set that describes the plane's destination. 
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FV This is set to .true, if the q file is in finite-volume format (i.e., the q variables are 
evaluated at the cell centers and not at the grid points). See section 2.4.2 for a more detailed 
description of the finite-volume option. 

GEOM This parameter should be used when a wall boundary is defined by high surface 
gradients and the standard grid redistribution has moved points in such a way that the 
original shape has been deformed. When GEOM is set to .true., the code will add the surface 
curvature function to the flow-field gradient function in the wall boundary regions. Points 
will thus be maintained in or redistributed into regions of high surface curvature as well as 
into regions of high flow-field gradients. The contribution of the geometry function will 
proportionally decrease away from both boundaries, with the internal lines controlled by the 
flow field only. If GEOM=.true. and QFUN=.false., adaption will be to geometry gradients only, 
for all grid lines. 

IJPLANE, IKPLANE, JKPLANE These parameters define the plane on which adaption 
takes place. By default, they also imply the direction of the stepping plane. The input 
parameter file needs only one of these parameters to be set to .true., and the code will 
automatically assign .false, to the other two. IJPLANE=.true. indicates that the plane 
represented by (i,j) will be the adaption plane and that plane stepping will occur in the k 
direction. Whether k is a forward or backward step will depend on KST and KEND (i.e., if KST 
> KEND then backward stepping will occur). Similarly, IKPLANE=.true. indicates that the 
adaption plane contains the (i,k) directions and that j is the plane-stepping direction. Finally, 
JKPLANE= . true . refers to the plane containing the points (j,k) and i is the plane stepping 
direction. 

IMPORT For multiple grids only. An input parameter list containing IMPORT=.true. 
must immediately follow a parameter set containing EXPORT=.true. The IMPORT list 
describes the destination (i.e. the receiving plane) of the plane of data defined in the 
EXPORT=.true. input-parameter set. The use of MGPLS is important for the subsequent 
adaption of the import grid. 

INDQ This parameter indicates which of eight possible flow-field variable(s) will drive 
the redistribution of grid points. From the standard PLOT3D format, five options are available: 

1 — » density, p 

2 — > x-momentum, pu 

3 — > y-momentum, pv 

4 — > z-momentum, pw 

5 -» stagnation energy, e 

Three more options are available by setting INDQ=6, 7 or 8 

6 -> pressure, p 

7 Mach number, M 

8 -» temperature ratio, T 

Pressure, Mach number, and temperature ratio are computed using the ideal gas 
relationship and assuming the Q file contains the standard variables. (The code actually 
assigns pressure to NDIM+1, Mach number to NDIM+2, and temperature to NDIM+3, so if the 
user has changed the value of NDIM to accommodate extra flow-field variables, INDQ must 
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reflect this change.) Note that INDQ=4 is not available for 2-D datasets. The user will normally 
choose to adapt to the flow-field variable that most represents the flow features. However, if 
different variables demonstrate different features, it may be advantageous to combine them to 
bring out all the features on the adapted grid. In this case, set INDQ=0 and input values for IQ. 

INTER indicates whether to use a linear (INTER=2), a quadratic Lagrange polynomial 
(INTER=3), or a cubic spline (INTER=4) scheme for interpolations. Interpolation is used 
throughout the code; for example, the q values in the output function file are interpolated at 
the new adapted grid points. Linear interpolation will usually provide the appropriate result. 

IQ is an array of eight (NDIM+3) integer values that are used only when INDQ=0. They 
allow the user to modify the adaption variable to a combination of variables. The order of IQ is 
consistent with the order of the flow-field variables (i.e., IQ(2) represents pu and IQ(7) 
represents M). The value of an index is the proportion that the corresponding variable will 
contribute to the final adaption variable. For example, IQ(1)=1,IQ(7)=1 [i.e., (1,0,0,0,0,0,1,0)] will 

produce an adaptive function of -(— + —) and IQ(1)=1, IQ(3)=1, IQ(6)=3 will produce 

2 ds 

+ + Obviously, (1,0,0,0,0,0,0,0) is the same as INDQ=L 

5 3s 5 dj 5 ds 

1ST, IEND contain the indices defining the first and last boundary lines of the adaptive 
domain in the i direction. Similarly, JST, JEND define the domain in the j direction and KST, 
KEND define the domain in the k direction. These variables define the limits of the adaption 
domain and must lie within the input grid boundaries defined for the physical domain,i.e., 

1 < IST,IEND < IMAX; 1< JSTJEND < JMAX; and 1 < KST,KEND < KMAX. 

Forward and backward stepping are also controlled by these parameters. If plane stepping is in 
the k direction, setting KST > KEND will produce backward stepping. Similarly, if stepping 
within the plane is in the j direction, setting JST > JEND will produce backward stepping. The 
reversing of the data is handled internally and is imperceptible to the user. Reversing either of 
the stepping directions will completely change the resulting grid, since the redistribution along 
the initial line and plane will be different, as will the connecting torsion springs. Reversing 
the order of the adaption direction (i.e., i in the default case) should have no effect on the 
solution, since the solution along a line is independent of the order of points. However, for 
meshes that contain very large and very small As,, numerical accuracy may be influenced by 
the adaption direction. 

ISTEP, JSTEP, KSTEP These are used in conjunction with the PLANE parameters 
described above and define the marching and adaption directions within the defined plane. 
Only one of these three parameters should be input and set to .true., and the code will assign 
.false, to the other two. If IJPLANE=.true., then only ISTEP or JSTEP can be true (KSTEP must 
be false). If ISTEP=.true., stepping occurs in the i direction and thus the adaption direction will 
be j. Points will be adapted along each constant i line and stepping will occur to the next i line, 
forward or backward, depending on 1ST and IEND. If JKPLANE=.true., only JSTEP or KSTEP 
can be true and similarly, for IKPLANE=.true., only ISTEP or KSTEP=.true. will have any 
meaning. The code puts out an error message if these inputs are inconsistent. 

IS, IE, JS, JE, KS, KE These parameters are used only for an export or import transfer 
process in a multiple-grid file. These variables define the range of the transfer domain within 
the plane, MPLANE. The PLANE parameter (e.g., IJPLANE) defines which coordinate plane is 
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being transferred. As an example, if IKPLANE— .true., then IS, IE, KS, and KE are used to define 
the domain within the transfer plane. If they are omitted, it is assumed that the entire plane is 
being transferred. 

LNSING In some 3-D grid types, planes emanate from a common grid line, generally 
on a wall surface. If the chosen set of adaption planes includes this common line, then this 
line should be adapted only on the first plane and not on subsequent planes. This adapted line 
is then placed in the first line of all planes. Although the input option is LNSING = «, it is 
likely that rt=l. 

LSTADD, LEND ADD These are input only if ADD * 0 and if only a portion of the grid 
is to be expanded. If they are not input and ADD * 0, the entire adaption domain (not physical 
domain) is assumed. If ADD=0, their value is ignored. 

LSTSUB, LENDSUB These parameters are input only if SUB * 0; they define the limits 
in which points are to be removed. If they are not input, then the entire adaption domain is 
assumed. 

MARCH This parameter refers to stepping within the plane. If the last line to be 
adapted ( j tnd ) is within the physical grid boundary (i.e., j eru i<j m!a )> a sharp discontinuity will 
occur between the last adapted line, j eru j, and the nonadapted line, j en d+\- Setting 
MARCH=.true. causes the remaining lines within the plane (i.e., j en d + i~* jam) to be realigned 
so that they are proportional to the last adapted line. This realignment will be performed for 
every plane. The process is not an adaption but has proved to be a useful tool. 

MARCHPL This parameter refers to the plane-stepping direction, and can be used 
independently or in conjunction with MARCH. If the last adaption plane is within the 
physical boundary (i.e., k^Kk^), each line in each subsequent plane will be proportioned 
with respect to the adapted lines in the k tnd plane. 

MG1, MG2 When NEDGE is non zero, an override mesh spacing is computed at the 
requested boundaries (either first, last, or both). This edge-point spacing, which is not a 
function of the adaption but of the initial grid, is merged into the three adjacent points, to 
produce a smooth transition. The default value of MG1 and MG2 is zero when NEDGE — 0, but 
four when NEDGE is requested. The user has the option of overriding this value, setting it to 
any other integer value. This can be used when adapting a boundary layer; when MG is 
increased, the dense edge spacing is maintained over a larger region. 

MGPLS This input variable is analogous to MGSTEPS described below, but applies to 
planes. If the first adaption plane (k st ) is internal to the physical domain, MGPLS = n permits 
the user to gradually bring in the effect of the adaption parameters to produce a smooth 
transition across the start plane. No adaption will take place on the first plane ( k sl ), and after n 
planes, full adaption occurs with the adaption parameters C',X, etc. equaling their input 
value. This is an especially important feature for grids that have an initial distribution on the 
wall boundary that defines a physical shape and cannot be changed. It is also useful for 
retaining matching multiple-grid boundaries, (see IMPORT) 

MGRID This parameter is only used for multiple-grid formatted files and indicates 
which grid to adapt. Do not use MGRID=1 for a single grid: setting MGRID to non-zero 
automatically implies a multiple-grid file. 


26 


MGSTEPS (m ? ) provides continuity when the first adaption line on each plane is 
internal to the physical boundary. Inputting MGSTEPS = n tells the code to start with no 
adaption on the initial line (i.e., retain the original distribution on the j s , line) and to linearly 
increase the adaption effect until, after n lines, full adaption occurs. At this point, C,, X, etc., 
will coincide with their input values. If MGSTEPS = 1, no adaption will be performed on the 
first line, but full adaption will occur on the second line. Figure 9 shows an example with 
m g = 5, and section 1.4.1 describes the parameter in detail. 

MPLANE Used only for multiple grids during an import or export process. MPLANE is 
an integer defining the transfer plane. The PLANE parameter is used (e.g., IJPLANE=.t.) to 
indicate the coordinate direction of MPLANE. 

NEDGE is a flag that requests an override on the computed side-edge boundary spacing. 
Side edges occur at i s , and i end and frequently need special handling. If there are no flow 
gradients near the edge of the domain, the standard adaption algorithm will pull points away 
from the edge. This may not be a satisfactory result, as, for example, when the side edge is 
internal to the physical grid boundary. Figure 10(a) shows a side-edge adaption with NEDGE=0. 
The flow-field gradients are concentrated in the center of the grid, and the first adaption point 
(i=4) has been pulled far from the boundary line at i sl = 3. It is clear that it is preferable for the 
adapted side-edge spacing to be continuous with the juxtaposed spacing in the external region. 
Even if the two boundaries coincide (i.e., i st = 1), the user may prefer a different spacing than 
that computed by the adaption algorithm. In either case, NEDGE can be set and the code will 
try to improve the side-edge spacing. The computed mesh-size override is merged into the 
next four points, but this number can be changed by MG1 and MG2. The result of setting 
NEDGE=1 is shown in Fig. 10(b). Depending on the case, both or only one of the side spacings 
may need improving. NEDGE=1 requests both edges, NEDGE=2 requests start edge only, and 
NEDGE=3 requests end edge only. 

NFILT is a "filtering" variable that sets the number of passes used to smooth the 
gradient of the input q data and the computed tension parameter, ©. The default value of two 
will generally suffice, but if the flow-field variables are discontinuous, it may be helpful to 
increase the value of NFILT. An increase in NFILT can also be used to expand or spread out a 
very sharp flow feature. 

NOUP This variable is used to increase or decrease the number of points in the grid (by 
using ADD or SUB) without performing an adaption. NOUP stands for NOUPdate. 

NOQ If no q file is available, SAGE can still be used to smooth the grid: perhaps to 
equal spacing or to the geometry function. If NOQ = .true., no file will be read on unit 8, and 
instead a constant flow field will be generated internally. This is quite different from QFUN = 
.false, where the q variables are not used, but are interpolated onto the new grid and output on 
unit 11. 

ORTHS, ORTHE The code assumes that orthogonality to the marching boundaries is 
desirable. This may not be the case, as, for example, in outgoing flow where the shape of the 
outer boundary is arbitrary. Setting ORTHS(l) = .false, will turn off orthogonality from the first 
to second adaption lines, and ORTHE(l) = .false, performs the same function when 
approaching the end adaption line. ORTHS(2) and ORTHE(2) will similarly affect the plane 
boundaries. 
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PLSING This is a parameter unique to 3-D grid adaption. It stands for plane singularity, 
and implies that a plane n is actually a single line, but is stored as a set of identical lines. The 
code will not be able to compute normals and will certainly "blow up" unless informed of this 
condition. It is only relevant when the adaption plane direction coincides with this collapsed 
plane. Do not use it when adapting in other plane directions. 

QFUN This parameter permits the user to adapt to the geometry function only. When 
QFUN = .false, and GEOM = .true, are input, the coefficient of the flow-field gradient (C q ) is set 
to zero for all grid lines (see section 1.5). In this case, the geometry function is computed 
throughout the grid and drives the adaption for all grid lines. 

RDSMAX, RDSMIN control the density of the redistributed points and are the 
maximum and minimum allowable grid spacings. They are input as proportioned values and 
are changed to physical variables internally, i.e., As max x s ^ /(n { ■ - \) and x j m(U . /(«, - 1), 

where n, is a constant equal to the total number of points along the adapted line, and is 
the length of the current adaption line. This implies that RDSMAX > 1.0 and RDSMIN < 1.0. (If 
both are set to 1.0, a uniform grid will result.) As an example, RDSMIN = 0.5 will prevent a 
converged As from being less than half the average step size. (Since many factors influence the 
distribution of grid points, this control is not absolute.) Note that since the adaption along the 
first line is not influenced by the torsion parameters, this initial line will present a clearer 
picture of the effect of RDSMAX and RDSMIN. 

REMOVE It is possible that an initial grid has unnecessary grid points in the outer 
region. Once an initial solution has been obtained, the user can see that these points are 
wasted. REMOVE = n will remove n points from the end of the adaption line, before 
performing the adaption. Use NOUP if only removing points is required. To remove points 
from the inner region, reverse 1ST and IEND. This will not change the adaption but will fool 
the remove operation. Remember, REMOVE deletes the boundary points, SUB retains the 
outmost one. 

SAVE This is useful when more than one adaption pass is made in the same program 
run, for example, an adaption stepping in the j direction followed by one stepping in the i 
direction. The output grid and function files are large, and setting SAVE=.false. on a set of 
$NAMEL will suppress the output for that particular adaption. If SAVE=.true. (default), each 
subsequent output set of xyz.out and q.out files will be assigned to different unit numbers. As 
stated in the execution section, the first output set is assigned to units 10 and 11. The second 
output set will therefore be units 12 and 13, and so on. Note: SAVE=.false. cannot be used for 
multiple grids. 

SUB When SUB = n (*0), points are removed from the adaption line. As an example, if 
n -l , every other point is deleted. If n=2, two consecutive points are deleted between the points 
retained. Note that the number of stepping lines remains constant: points are only removed 
from the adaption line. To remove points from both directions, two passes are required. See 
LSTSUB, LENDSUB if only a selected range of deleted points is to be deleted. 

TWOD If the input grid and function files are stored as 2-D PLOT3D files, this 
parameter must be set to .true. The code will assume that IJPLANE is the adaption plane and 
JSTEP the stepping direction. The user must specify ISTEP=.true. if required. The next section 
discusses the adaption of 2-D problems. 
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2.4 Alternative Grid Types 

2.4.1. Two-dimensional Adaption. The SAGE code can accommodate 2-D datasets, 
and will adapt the single plane in the same manner as in the original 2-D SAGE code (Davies 
and Venkatapathy 1989). When the input parameter TWOD is set to .true., the code will read 
the grid (assuming ( x,y )) and function files as 2-D files. These datasets will then be internally 
converted to a 3-D format; the number of k planes will be set equal to one, creating a constant z 
coordinate, and the 4th q function will be shifted to index 5. (Note that indices INDQ and IQ 
retain their 3-D relationship). Because of this reorganization, no special handling of 2-D 
datasets is required within the body of the code. Datasets are reconverted to 2-D form before 
output. To accommodate larger dimensions, the parameter statement may be changed at the 
beginning of each routine. Since 2-D datasets require only KD=1, ID, JD and IMX may be 
significantly increased. However, this change is not made automatically and the user must 
change the parameter statements if necessary. 

2.4.2. Finite-volume Grids. The solution file associated with finite-volume 
applications contains q values evaluated at the cell center. Therefore the IMAX etc. values on 
the header record are one less than the values given in the grid file. If SAGE finds that the size 
records disagree and that the finite volume option is off, an error message is sent to the user. 

Care must be taken with the boundary (or ghost) cells. SAGE interpolates q onto the 
internal grid points and then sets all boundary values (in the physical domain) equal to the 
adjacent interior value. After adaption, flow values are interpolated back to the cell center. If 
all planes are not adapted, the q values in the final adapted plane will be interpolated as if they 
are on a 2-D surface since there will be a discontinuity between the adapted and non-adapted 
grid points. The boundary values of the physical domain are again set equal to the adjacent 
values, regardless of whether the entire physical domain has been adapted. It is therefore very 
important for the user to check all the boundary cells in each coordinate direction. 

Related to the handling of the boundary cells, SAGE will only adapt a finite-volume 
grid if it is a single surface (whether defined in 2-D or 3-D) or is a 3-D grid with four or more 
planes. 

2.4.3. Multiple Grids. The multiple-grid format is a single file containing a 
collection of separated grids. Preceding the grids are two header records, one describing the 
number of grids and the other the size of each grid. The associated q file is similarly defined. 

Some complex multiple grids utilize the blanking feature available in PLOT3D. At this 
time, adaptions by the SAGE code cannot take into account these blanked regions. If 
overlapping regions are adapted in separate grids, it is the user’s responsibility to interpolate 
the results. However, matching zonal planes can be handled with the plane-transfer feature 
described below. 

2.4.3.1 File handling. Since each grid is stored sequentially within a multiple-grid 
file, SAGE reads and copies all grids (and their associated solutions) to the appropriate output 
files until the requested grid (MGRID=n) has been read. This grid is now adapted (or a plane 
transferred) and written, along with the interpolated flow solution, into the correct sequence 
in the output files. Finally, any grids following the adapted grid are also read and copied to the 
output files. Subsequent adaption or export/import passes will use these output files as input 
files. If several passes occur in the same computer run, several sets of large files could be 
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created and the user is reminded to delete unneeded files. The SAVE parameter cannot be 
used. 

An additional set of files will be created if any grid has changed size (e.g., the ADD or 
SUB option has been used). Grid sizes are stored in the header record that has already been 
written to the output file before adaption took place. This header record must be amended to 
reflect the new grid size. The final output files are therefore read in as input files, the header 
record is corrected, and yet another set of output files are created. An output message keeps the 
user informed of the unit numbers for the final set of output files. 

2A.3.2. Data transfer using EXPORT and IMPORT. Section 1.7 explains the need 
for data transfer between grids. A set of input parameters are available that control the transfer 
of data from one domain within a specified plane of a specified grid to a matching domain 
within another grid. The same namelist format is used that controls the normal adaption 
procedure in SAGE, but two sets are required in the user-input file: the first describes the 
'export' plane and the second describes the 'import' plane. Here, the term ’export' means a 
domain (i.e. a surface, plane or sub-plane) whose data will be transferred. The domain 
receiving this data is the import' domain. A transfer domain is defined by 

1. the transfer code, either export-.t. or import=.t. 

2. the grid number, mgrid 

3. the plane direction, ( ijplane,ikplane,jkplane ) 

4. the plane number, mplane (this will often be the first or last plane) 

5. the range of points describing the domain (i.e. the start and end points in two 

directions) within the plane 

Notes: (a) For 2-D datasets only one line in a plane is transferred and ijplane =.t. 

and mplane=l are defaulted by the code. The transfer domain is described by is, ic, js, je and 
one of these directions must have equal start and end points. If the 'export card is the first 

input set, remember to include twod=.t. 

(b) Multiple transfers can be handled in the same run of the code as well as 

a combination of adaptions and transfers. See the examples section for clarification. 

2.5 Output Message File 

The sage. out file is written to unit 6, the normal default for the output screen. It 
contains messages that help explain what has happened during program execution. At the end 
of each adaption, "ADAPTION n COMPLETE" indicates that the program was able to run to 
completion and that xy.out and q.out files have been created. The message "OUTPUT FILES 
ON UNITS n, AND n 2 " informs the user where the final set of files can be found. 

The following are other messages that may be seen (given in alphabetical order) along 
with a short description of their meaning. 

ADD OPTION EXCEEDS DIMENSION (Critical) 

Self-explanatory. Increase array dimensions. 

FINITE VOLUME METHOD NEEDS 1 OR 4+ PLANES (Critical) 

A single plane is treated like a 2D surface: i.e., 4 points are used for cell-centering. Four 
planes are needed for the finite volume method in 3D. 
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GRID HAS IDENTICAL POINTS AT i AND i+1 ON LINE; AND PLANE k : 

CHANGE ADAPTION DOMAIN OR REMOVE SINGULARITIES (Critical) 

Computation of body normals is impossible in this region. 

GRID SIZE TOO LARGE FOR DIMENSION STATEMENT: 

CHANGE PARAMETER STATEMENTS TO lD-a,,JD=n 2 , KD=n,, IMX=n,: 

THESE ARE THE MINIMUM DIMENSIONS FOR THIS ADAPTION PASS ONLY (Critical) 
Increase array dimensions to size suggested. Note that these values are appropriate for this 
set of input-control parameters only and may need to be changed for other adaptions of the 
same grid. For multiple grids, this message may be preceded by additional information. 

MPORT CARD EXPECTED, NOT FOUND (Critical) 

The previous card in the input stream contained export=.t. This must be followed by a set 
containing import=.t. 

INCONSISTENT PLANE AND STEP (Critical) 

A stepping direction has been requested that is not available for the requested plane. For 
example, IKPLANE and JSTEP. 

INPUT FILE SIZES DO NOT MATCH (Critical) 

The grid dimensions on the header records of the grid file and solution file do not match. 

IS THIS FINITE VOLUME? IF SO, SET FV=.TRUE. (Critical) 

The grid file dimensions are one greater than the solution file. Should this be finite 
volume? 

MAX I TOO LARGE, CHANGE ID TO n, 

MAX J TOO LARGE, CHANGE JD TO n 2 
MAX K TOO LARGE, CHANGE KD TO n } 

CHANGE IMX TO n 4 (All critical) 

These messages occur only for multiple grids. It implies that one or more of the grids in 
the file is too large for the dimension statement. It need not be the adaption grid. It is 
possible that even if this is corrected, a subsequent run could indicate a second dimension 
change to accommodate data swapping for the adaption grid. 

NDIM TOO SMALL FOR 2D TO 3D TRANSFORMATION (Critical) 

Increase NDIM dimension (e.g., from 4 to 5) so internal transformation can be made. 

NO CONVERGENCE ALONG INITIAL LINE, ERRMIN=a, (Warning) 

The initial line is a 1-D adaption only. This is rarely a catastrophic error, especially if a, is 
small; however, the adaption may not be completely satisfactory. The only control 
parameters that affect the initial line are RDSMAX, RDSMIN, and NEDGE. 

NO CONVERGENCE ON LINE ; AND PLANE k, ERR=a 2 (Warning) 

This message is only a warning and adaption continues. Even many of these messages 
may of be no concern as long as a 2 is small. If adaption is successfully completed, check the 
new mesh to see if it is acceptable. 

NO POINTS ADDED (Warning) 

Inconsistency in parameters requesting adding points. 

NO OUTPUT FILES (Warning) 

For whatever reason, no files have been output. (Is save=.f. on all passes?) 
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NUMBER OF POINTS INCREASED FROM n, TO n 2 (Informational) 

If the ADD option has been input, this message gives the new grid dimensions. 

NUMBER OF POINTS DECREASED FROM n, TO n 2 (Informational) 

If the SUB option has been input, this message gives the new grid dimensions. 

PLANE n,, GRID n 2 COPIED TO PLANE n^GRID n 4 (Informational) 

The successful transfer of a surface from one grid to another in a multiple grid file. 

POINT(S) REMOVED, NUMBER OF POINTS NOW n (Informational) 

REMOVE option has been invoked. 

S IS NOT MONOTONIC ON LINE j AND PLANE k (Critical) 

This message will terminate the program. It indicates that the values of s, at the 
completion of the iteration on line j are not monotonically increasing, thus implying 
crossover of points. Since this is unacceptable, the program outputs the data. It is then 
possible view the plots and reevaluate the control parameters. 

SOLUTION FILE IS F-V: CHECK YOUR BOUNDARY CELLS! (Informational) 

Data in ghost or boundary cells may be incorrect. 

SUB OPTION PRODUCES TOO FEW POINTS FOR ADAPTION (Informational) 

Fewer than 10 points remain in the adaption direction. 

TOO FEW POINTS FOR ADAPTION WITH NEDGE=1 (Critical) 

There are less than 10 points along the adaption line. This is not appropriate, especially if 
nedge is set. 

WARNING: DIRECTIONAL SIGNS INDICATE COORDINATE SYSTEM 

MAY NOT BE RIGHT-HANDED (Critical or informational) 

The calculation of body normals assumes a right-handed coordinate 

system (an example is shown here, where j goes into the paper). In most 
cases, this message indicates the grid must be re-ordered. 

2.6 Outline of Each Subroutine 

The MAIN routine is a driver routine whose task is to call the relevant subroutines. A 
loop (using NADS) is set to provide for multiple adaption passes. The first routine to be called 
is INITIAL, which reads and organizes the data. A loop is then set up for the adaption of each 
plane, and the constant planar variables are computed. Finally, a loop is set up for each line in 
the plane in which all the coefficients of the adaption equation are computed, followed by the 
solution process. When each pass is complete, the OUTPUT routine is called. 

Along with the MAIN routine, the SAGE code consists of the following subroutines, 
listed here (as they are in the source code listing) in alphabetical order. Arguments shown in 
boldface are computed within the subroutine. 

ADDPTS( IER) 

This routine is called by MAIN if ADD * 0 . Extra points (depending on ADD) are inserted 
between every grid point within the range LSTADD to LENDADD, by a linear interpolation. 
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.ADDV(COSXl,COSYl,COSZl,Al,COSX2,COSY2,COSZ2,A2,COSX,COSY,COSZ) 

This is a utility routine. The direction cosines of two unit vectors (COSXl,COSYl,COSZl) and 
(COSX2,COSY2,COSZ2) are input arguments. The routine computes the direction cosines 
(COSX,COSY,COSZ) of the unit vector that represents the sum of the input vectors, 
proportioned by the coefficients A1 and A2. 

BLOCK{ J,K) 

Initially, all grid coordinates are stored in the input X, Y, and Z arrays. In this routine, a block 
of data around the current j line is stored in arrays XJ, YJ, ZJ to prevent the original grid from 
being overwritten during interim calculations. Only the converged adaption line is replaced 
into the original grid arrays. The XJ, YJ, and ZJ are made up of all i points on the current j line 
and all i points on the j-1 and j+1 lines in the k-1, k and k+1 planes, giving (n f , 3,3) 
dimensioned arrays. At boundaries, nonexistent data points are filled with 999. Most 
calculations within the code, but especially the computation of the normal vectors, are 
performed on this block. 

CROSS V(XT,YT,ZT,XT1,YT1,ZT1,DST,COSV,AAP,DAP,ICROSS,J) 

Appendix II, section 1.9.2 contains the analysis used to develop this routine. COSV is the array 
containing the direction cosines of the vector v , defined in that appendix. (XT,YT,ZT) are the 
coordinates of a j line, (XT1,YT1,ZT1) are the coordinates of a juxtaposed line, and DST is As 
along j. As an example, this routine is used to compute s - s' (i.e., AA') and DA' , the distance 
between s' and the corresponding node at (i,j-l,k), as shown in Fig. 4. ICROSS is the array 
containing / and indicates the intersecting segment for each COSV. 

CSPLIN( NT,S,V,SPF) 

The cubic spline coefficients, SPF, are computed for NT points. S is the streamwise location of 
the function V. These coefficients are used by the routine SPEVAL if GEOM=.true. or 
INTER=4. 

DETERM( A1 ,B1 ,C1 , A2,B2,C2, A3,B3,C3,DET) 

This routine computes the determinant of the three vectors whose direction cosines are given 
in the argument list. 

DLENGQ L,K) 

When NEDGE is set, the tension parameter co is amended at the edges (i.e., at 1ST and IEND) to 
improve edge-boundary spacing. This routine computes DLENGS and DLENGE, which are 
used in the edge co calculation (by the routine WTEDGE). The values of DLENGS and DLENGE 
depend on whether the grid is defined outside of the adaption domain. JL indicates which j 
line is needed. 

EDGEMGCV AR) 

This is a utility routine. For various reasons, the values of some variables at the 1ST and/or 
IEND edges are overridden. To blend these different values into the calculations, this routine 
will merge the new values (given at the two boundaries of VAR) into the next three (or MG1, 
MG2) grid locations of VAR. 
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FBAR(J,K) 

This routine is called once for every j line. The As and the gradients dq/ds are computed at the 
input grid points and stored in FQ. If duplicate points are found (i.e., 3 a As, = 0) a message is 
printed and the program terminates. If GEOM=.true., the wall gradients dg/ds are also 
computed by calling WALLS, and stored in FG. In addition, the coefficient C g is computed and 
FG is added to FQ and stored in F. The routine INTF is called to interpolate the value of F at 
these new grid points and to compute the normalized form of F, i.e., FB. The routine GETB is 
called to find the value of B for this j line. For lines other than the first, initial guesses for the s 
distribution are extrapolated from the converged j, along the j-1 line. Since these do not 
correspond to the input points, new local values of x, y, and z are interpolated by calling 
PROPS. 

FILTERS AR,NIPTS,NFILT) 

This routine smoothes the parameter contained in VAR by adding a second derivative term, 
v, =.75v,+.125(v l+1 -v f _,h NFILT is the number of smoothing passes (default=2). The parameters 
smoothed are /,■ = f (dq/ds) and to,; they are returned to the calling routine via VAR. 

FVORG(IND) 

If IND=1, this routine interpolates the cell-centered q values onto the grid points to enable 
adaption to take place in the normal manner. Values at the first and second boundary points 
are set equal. If IND=2, the q values are reinterpolated to the new cell centers. 

GETB(J,K) 

This routine computes the value of B used to evaluate os. B is found by an iterative process and 
is said to converge when the minimum requested As equals the computed minimum As . The 
analysis for this routine is given in appendix I. 

GETWT 

This routine computes co ( , the modifier of to that is set when any computed As lies outside the 
requested range. 

HEADIOaC, IER) 

This routine reads (IC=1) and writes (IC=2) the header records for multiple-grid files. It also 
tests to see if the grid dimensions are too large for the programmed dimensions. 

fN/TML(NOMORE,MTCH) 

This routine sets all the default values and reads the input parameter file. If EXPORT is true 
for a multiple-grid file, MTCH is set to .true, and control is handed to the MATCH routine. In 
all other cases, the appropriate grid and function file input routines are called (either READAT 
or READMULT). When necessary, the grid points and corresponding flow-field data are 
rearranged to correspond to the data organization assumed by the analysis. This routine also 
controls the addition or removal of points. If no more adaptions are requested, NOMORE is set 
and control is returned to the main routine. 

INTF( F1,F2,S1,SMID,NPTS) 

This routine interpolates to find F2 (new F) at the new s, (SI), based on the input values of FI 
(current F) and the midpoints (SMID) of the input s array. The interpolation routine chosen is 
based on the value of INTER. 
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/NTXYZQ(J,K,J1,K1,SS,SN, QJ) 

This routine interpolates for X, Y, Z, and Q at the new SN, given the corresponding values at 
SS. The new (x,y,z) coordinates are stored in the block of data defined by XJ, YJ, and ZJ. The q 
data is stored in QJ. The appropriate interpolation routine is called, its choice based on INTER. 

LAGCOF(SNEW,SARR,NPTS,M,Pl,P2,P3) 

This routine computes the Lagrange coefficients PI, P2, and P3 for a point SNEW, with respect 
to the input s array, SARR. These are used by the calling routine to interpolate for the variable 
at SNEW. First or second order is available; the choice depends on INTER. M is the associated 
index for PI and will reflect a forward or backward interpolation, depending on the location of 
SNEW within the interval. 

LINE1 

This routine solves for the adapted values of s, along the initial line j = j sl on the initial plane 
k s[ . Since both torsion terms are zero on this line, the As, are computed from the 1-D approach. 

MARCHJ( K) 

If the final adapted line within plane k is internal to the physical end-boundary line, this 
routine will redistribute the points on the remaining j lines, based on the distribution along 
the j end line. The routine is called for each plane. This action is performed only on request by 
the user and is not an adaption to the flow field. However, it will prevent the discontinuity 
between the last adapted line and the remaining non-adapted grid lines. 

MARCHK 

This function of this routine is similar to MARCHJ but it is called on user request, if the final 
adaption plane k end is internal to the physical end plane k^. It will redistribute points on all 
lines on the remaining planes to be proportional to the corresponding line on the last adapted 
plane. 

MATCH(NOMORE,IER) 

This routine is called as soon as an 'export' card is read. If necessary, the output files are 
rewound ( REWND ); and then the header records are read ( HEADIO ). The export grid is read 
into core ( MULTIO ) and the requested plane stored ( STOREX ). Now, the next set of input 
parameters are read (which should be an 'import' set). The input files are closed, rewound and 
re-read up to the import grid ( READMULT ). The stored export plane is copied into the import 
plane ( STORIM ) and all the files are written to the output units (WRIT MULT). 

MGWALLS( J) 

This routine is called when GEOM=.true. It computes the coefficient of the geometry function, 
C g (FGW), based on the local aspect ratio. 

MULT/0(IN1 ,IN2,IOUTl ,IOUT2) 

This routine reads and immediately writes multiple grid files. INI and IN2 give the range of 
grids to read in (e.g., if MGRID=3, READMULT will initially call MULTIO with IN1=1 and 
IN2=3). IOUT1 and IOUT2 give the range of grids to output (in the same example, 
READMULT will set IOUTl=l and IOUT2=2, since grid 3 will be written after adaption). The 
header records are handled separately in HEADIO. 
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NOADAPT(],K) 

This routine updates variables and/or places them in appropriate arrays when no adaption is 
to be performed on the current j line, but adaption is to be performed on the next line. In this 
scenario, the code expects certain variables to be available at j-1. This occurs, for example, in 
merging situations (MGSTEPS > 0) and for common lines (LNSING > 0). 

NORM(Fl,NPTS) 

The function FI is normalized as (FI - Fl^KFl^ - Fl mm ). If maximum and minimum 
values are equal or nearly equal, the normalized variable is set to 0( 10 -5 ). 
NORMPT(IP,JP,KP,INDPL,PA,PB,PC,SING) 

This routine finds the vector at the point (IP,JP,KP) normal to the plane defined by INDPL. 
Although three direction planes exist through the point, only two are needed by the calling 
routines: INDPL=1 is the (i,k) plane and INDPL=2 is the (i,j) plane. The analysis to describe this 
routine is given in Appendix III of this report. Four normals are computed and the average is 
found, with the direction cosines returned in (PA,PB,PC). SING is set to 1 if the normal does 
not exist. 

OinTUT(NADS) 

OUTPUT is called at the conclusion of each adaption set. Proportioning of any remaining 
planes is performed if requested by MARCHK and then the data files are returned to their 
original order, to conform with the input mesh structure. Either WRITMULT or WRITOUT is 
called to output the grid and flow-field files. NADS indicates which adaption number has been 
completed (maximum of 10 sets). This routine is also called if s becomes nonmonotonic 
(OK=.false.). In this case, the new mesh points that have been computed are output to help the 
user choose more appropriate control parameters. 

PROPSQ, K) 

After a new solution of s, has been obtained on a line /, the code stores this data in line j-1 and 
steps to the next line. This routine proportions the new j, throughout the non-adapted 
regions of the block of data defined by XJ,YJ,ZJ. Essentially, this provides a first guess for the 
current j line and also produces smoother planes for the vector normal calculations. (See 
Appendix III) 

Pl/RPLE( A1 , A2, A3,B 1 ,B2,B3, V 1, V2, V3,NFL) 

A and B are direction cosines of two vectors defining an enclosed plane. This routine takes 
their cross-product and normalizes the result to give the direction cosines (V1,V2,V3) of the 
unit normal to the enclosed plane. NFL is the direction sign of the normal. 

READAT 

This routine reads in the single grid and function files in PLOT3D binary format. Input 
datasets may be in 2-D or 3-D form. The size of IMAX, JMAX and KM AX on the header record 
is checked ( SIZING ) before the data is read. If NOQ is true, the Q array is filled with 1.0. 

READMULT(1ER) 

This is the read routine for multiple-grid files. Since some writing of files also occurs, the 
output unit numbers are verified. The header records are read and copied to the output files 
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( HEADIO ). Then MULTIO is called to read the requested grid into core and to write all 
preceding grid data onto the output files. Finally, if NOQ is true, the Q array is filled with 1.0. 

REWND 

This routine rewinds the output files and opens them as input files. 

SETUPR J,K) 

SETUPJ computes the direction cosines of the vectors u, b, and e used to evaluate the torsion 
vector F. These vectors are associated with the (i,j) points within the constant k plane. 

SETl/PK(J,K) 

This routine is similar to SETUP], but the direction cosines of u* , b *, and e * are computed to 

evaluate the plane torsion vector, F*. These vectors are associated with the (i,k) points within 
the constant j plane. 

SINGPLN 

When PLSING is set, SINGPLN stores the result of the first adapted line into every line of the 
same plane. 

SJZJNG(IER) 

The parameter statement at the beginning of each subroutine presets the dimensions 
(ID,JD,KD) of the grid and q files. This routine compares the input grid dimensions 
(IMAX,JMAX,KMAX) or (IM(NGRID),JM(NGRID,KM(NGRID)} to these preset values. If 
insufficient space has been allocated, the minimum possible values of ID, JD, and KD are 
computed for this adaption to proceed. (These values are not obvious since space must be 
allocated for data swapping.) IMX, the maximum of ID, JD and KD is also evaluated. A message 
is then sent to the user to recompile SAGE with the suggested dimensions. Finally, IER is set to 
1 to inform the INITIAL routine to terminate the code. 

SOLUTQ,K) 

By the time SOLUT is called, all variables have been computed that are needed to obtain the 
new distribution of s , . The coefficients of s t (see Eq. 8 in the first section of this report) are set 
up in a tri-diagonal matrix and solved for s l . Interpolated values of ©, are found at these new 
values of s, and iterations are performed until y is small or too many iterations 
have been performed. In both cases, a check is made to see if all the s, are monotonic. If so, the 
program continues; if not, the flag OK is set to false, causing the program to output the current 
grid and terminate. 

SPEVAL(NT,S,V,SPF,SI,VI,VPI,VPPI) 

This is a cubic spline interpolation routine used when INTER=4 or if GEOM is true. It uses the 
coefficients (SPF) computed in CSPLIN. In addition to interpolating for the variable V at SI, the 
corresponding first (VPI) and second (VPPI) derivatives are also evaluated. 

STOREX 

This routine stores the data from the export plane into XP,YP,ZP and QP. 

STORIM 

The export data in XP,YP ZP and QP are copied to the import plane. 
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SUBPTS 

When SUB = n±0,n points are deleted between each retained mesh point. If LSTSUB and/or 
LENDSUB are nonzero, the range of the deletion is restricted. Deletion occurs only in the 
adaption direction and does not decrease the number of stepping lines or planes. 

SWAP1NV 

This routine is called if k st > k end , j„ > j tnd , or i s , > i end . The order of (i,j,k) in the x, y and z input 
matrices are reorganized to ensure that internal computations have monotonically increasing 
indices. The flag for the handedness of the coordinate system is amended accordingly. 

SWAPXYZ( RSWAP) 

Since the internal computation assumes that j is the stepping direction within the plane and 
that it is the stepping direction of planes, this routine is called to interchange x, y, z, and q 
when the input requests alternative stepping directions. RSWAP is a flag that states whether 
this data exchange is at the start of the computations or is the reverse process required for 
output. This routine is lengthy due to minimizing storage requirements. 

SWAP2D(IO) 

TWOD=.true. indicates that the input datasets are formatted in two dimensions. This routine 
reorganizes the 2-D plane to appear in the code as a 3-D surface. Every z is given the value of 
zero, each Q(4) is moved to Q(5), and Q(4) is zeroed. IO indicates whether the translation is 
from 2-D to 3-D, which occurs on input, or from 3-D back to 2-D for output. 

rORCOF(L,JK,JKST,JKEND,MGNOS,MARCHJK) 

TORCOF is called twice in each loop, once with all the arguments representing the k plane 
passing through (i,j,k) and once with the arguments representing the j plane also passing 

through (i,),k)- This routine amends the coefficients (C,,X) etc. of the torsion vectors t and f 
based on the current line location. For example, C, is decreased when leaving or approaching a 
boundary to emphasize orthogonality. 

TORSIONQ.K) 

This routine first chooses the appropriate b and adds to u to obtain n. The torsion vector t is 
then obtained by adding it and e. V is computed in a similar manner. The routine CROSSV is 
then called to find the intersection of the torsion vectors with the j line from which both s' 
and s* can now be evaluated. Finally, a check is made to ensure that .v, and s ( monotonically 
increase. If they do not, the code attempts a correction, but any major problems will cause the 
code to terminate in the SOLUT routine. 

UM7V(X1,Y1,Z1,X2,Y2,Z2,DIRCX,DIRCY,DIRCZ) 

This is a utility routine that finds the unit vector from (X1,Y1,Z1) to (X2,Y2,Z2). DIRCX, DIRCY 
and DIRCZ are the direction cosines of this vector. 

UPDATED, K) 

This is the last routine called in the iteration loop for the current j. Newly adapted values of s, 
have been found. The values of x, y, z and q at this new distribution are interpolated 
(INTXYZQ) and replaced into the matrices containing the physical mesh. 
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V MERGE (DIRV,LST,LEND) 

This routine performs the same function as EDGEMG, but with a vector quantity (DIRV) in 
place of a scalar value. LST and LEND indicate which value of DIRV must be merged into the 
next three (or MG1, MG2) points. 

WALLS{ JW,K) 

Along line j (JW), the geometry gradient, as defined by the radius of curvature, is computed for 
each grid segment. Normally, j will equal j sl or j rnd . However, for cases when geometry is the 
only adaption variable (QFUN=. false.), WALLS is called for every j line. 

WRITMULT 

This is the write routine for multiple-grid files. The current adapted grid is now output and 
any subsequent grids are read and written by calling MULTIO. If the grid size has changed, the 
header record is updated; the output files are closed and opened as input files; the header 
record is corrected and all data are copied to the new output files. 

WRITOUT 

This routine writes the single-grid 2-D or 3-D adapted grid and interpolated function file on 
units NITG and NITQ. If the finite-volume option is set, IMAX, JMAX and KMAX are one less 
on the q file. 

WTEDGE( J,K) 

This routine is called by MAIN when NEDGE modification is requested. The edge values of 
the tension parameter at the next j line are a function of the average coA s along the just- 
completed adapted line. This routine calls DLENG to obtain the appropriate value of edge As 
on the next line, computes the average u>As on this line, and evaluates WDS and WDE to be 
used in routines LINE1 and SOLUT. 


2.7 Nomenclature 

The following is a list of variables used in the SAGE analysis. When applicable, the 
corresponding FORTRAN name used in the code is shown in boldface. 


A,B 

C * 

C, 

c; 

b' 

c c* 

d; ( d x , d y , d z ) 
e; (e x ,e y ,e z ) 


constants used to compute co, (A,B) 
aspect ratio, used to control C g , FGASP 

coefficient of fig) in co calculation, FGW 

coefficient of f(q) in co calculation, FQW 

input proportion coefficient for torsion, CT(1) 

input proportion for torsion between planes, CT(2) 

normal vector to j line within k plane; direction cosines of b , COSB 

normal vector to j line within j plane, COSBK 
modified values of C„ C,*, CTM(l), CTM(2) 

straightness vector, COSD 

average straightness vector within k plane, COSE 
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average straightness vector within j plane, from k 2 —* k l,COSEK 
gradient of q (and g if necessary), F 

minimum and maximum /used to normalize/, FMIN,FM AX 

normalized function of /, FB 
geometry function 

subscript indicating the current node in adaption direction, I 
total number of points in i, j and k directions of input grid file, 
IMAX, JMAX, KMAX 

start of adaption domain in i, j and. k directions, 1ST, JST, KST 

indices indicating end of adaption domain, IEND, JEND, KEND 

subscript of the current stepping line, J 

subscript of the current adaption plane, K 

torsion-related constant 

local subscript relating to node i, L 

number of stepping lines before full adaption, MG STEPS 

the plane equivalent of m g , MGPLS 

flag for edge control, NEDGE 

number of points in the adaption line, NIPTS 

merging coefficient for lines f(m g ), CNM(l) 

merging coefficient for planes /( m* ), CNM(2) 

orthogonality vector within plane, COSN 

orthogonality vector between planes, also stored in COSN 

input flow-field variable (p,pw,pv,pw,e), Q 

radius of curvature for geometry function, FGS 

stream wise length, SS or SN 

vector representing s, SSX, SSY, SSZ 

maximum value of s on line SMAX 

value of streamwise length used for torsion within planes, SP 
value of streamwise length used for torsion between planes, SPP 

requested minimum and maximum grid spacings, DSMIN,DSMAX 
computed minimum and maximum grid spacings 
torsion force 

within-plane torsion vector, COST 
torsion vector between planes, also COST 
proportion of b and u used to compute n, TN(1) 

proportion of b * and S* used to compute n" , TN(2) 
vector normal to j-1 line in k plane, COSU 

vector normal to j line on k-1 plane, COSUK 
input grid mesh, (X,Y,Z) globally and (XJ,YJ,ZJ) locally 
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X input magnitude of torsion control parameter, CLAM(l) 

X* input magnitude of plane torsion control parameter, CLAM (2) 

0 angle for torsion computation 

x within-plane torsion-related parameter, TAU 

to tension spring force, WEIGHT 

to, computed weighting on co, WT 

y between-plane torsion-related parameter, TAUPL 

2.8 List of Major Variables 

This section contains the list of variables (in alphabetical order) used in the SAGE code. Local 
variables that contain only intermediary values are not listed. The format is: 

Variable name(dimension) /r, /r 2 /brief description 

where r { describes what type of variable — input, local, parameter, or common block name 
and r 2 lists routine(s) where the variable is initialized. 


A 

AA(IMX) 

AAP(IMX) 

ACT 

ADD 

ALPHA 

AMACH(IMX) 


/ COM2 / INITIAL / A used to compute co 
/local/SOLUT/ coefficient of in solution matrix 
/argument/CROSSV/ s- s' or s- s' 
/local/TORSION/ final modified C, 
/COM19/input/ if set, add grid points 
/COM12/input/ information only for PLOT3D 
/ local /FBAR/ computed Mach number 


B 

BB(IMX) 

BCONV 

BJ1 

BK1 


/COM2/GETB/ B used to compute to 
/local/SOLUT/ coefficient of s, in solution matrix 
/ local /GETB/ convergence criteria for B iteration 
/local/GETB/ value of B along j-1 line in the k plane 
/COM15/GETB/ value of B along j line in the k-1 plane 


CC(IMX) 

CLAM(2) 

CLAMW(2) 

CNM(2) 

CONV 

COSB(IMX,3) 

COSBK(IMX,3) 

COSD(IMX,3) 

COSE(IMX,3) 

COSEK(IMX,3) 

COST(IMX,3) 


/local/SOLUT/ coefficient of s, +1 in solution matrix 

/COMlO/input/ X and X*, magnitude of torsion terms 
/COMIO/TORCOF/ modified CLAM 

/COMIO/TORCOF/ X modifiers for MGSTEPS* 0 and MGPLS* 0 
/COM15/INITIAL/ general convergence criteria 

/COM7/TORSION/ direction cosines of b , in (i,j) plane 
/COM20/TORSION/ direction cosines of b' , in (i,k) plane 
/local/SETUPJ,SETUPK/ temporary straightness vector, d 
/COM7/SETUPJ/ straightness vector, e , in (i,j) plane 
/COM20/SETUPK/ straightness vector, e* , in (i,k) plane 
/local/TORSION/ torsion vector, t or ?* 
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COSU(IMX,3) /COM7/SETUPJ/ u, normal to j-1 line, in (i,j) plane 

COSUK(IMX,3) /COM20/SETUPK/ u, normal to j line in k-1 plane 

C7(2) /COMIO/ input/ C, and C,*, directions for torsion vectors 

CTM(2) /COMIO/TORCOF/ modified C,and C t 


DAP(IMX) 

DAPPL(IMX) 

DET 

DLENGE 

DLENGS 

DMINSDB 

DS(IMX) 

DSMAX 

DSMIN 


/COM9/CROSSV/ DA' for s' calculation 
/COM9/CROSSV/ DA' for s’ calculation 
/argument /DETERM/ value of three-order determinant 
/COM6/DLENG/ As computed for end-edge modification 
/COM6/DLENG/ As computed for start-edge modification 
/ local /GETB/ dmin(As)/dB 
/COM3/FBAR,LINEl,SOLUT / As 
/COM6/FBAR/ As^ from RDSMAX 
/COM6/FBAR/ As m ,„ from RDSMAX 


EXPORT 


/COM21 /input/ if true, store data for transferring between multi-grids 


F(IMX) 

FB(IMX) 

FF(IMX) 

FG(IMX) 

FGASP 

FGS(IMX) 

FGW 

FQ(IMX) 

FQW 

FSMACH 

FV 


/COM2/FBAR/ flow gradient, / = dq/ds, at input s 
/COM2/INTF/ /, normalized / at current s 
/ local /SOLUT/ coefficient of RHS of solution matrix 
/COM17/FBAR/ normalized dg/ds 
/local/MGWALLS/ function of aspect ratio 
/COM17/WALLS/ dg/ds at SMSS 

/COM17/MGWALLS/ coefficient of FG for to calculation 
/COM17/FBAR/ normalized dq/ds 

/COM1 7 / MGW ALLS / coefficient of FQ for a> calculation 
/COM12/input/ information only, for PLOT3D 
/COM15/input/ flag to indicate finite-volume q file 


GEOM 


/COM11 /input/ request for geometry function 


ICROSS 

ID 

IEND 

IINVERSE 

IJPLANE 

IKPLANE 

IM(10) 

IMAX 

IMAXQ 

IMPORT 

IMQ(IO) 

IMX 


/ argument /CROSSV / index for interval location of s or s 
/dimension/parameter statement/ currently set at 75 
/ COM4 /input/ last node along i in adaption domain 
/COM13 /INITIAL/ adaption requested in backward i steps 
/COM14/input/ true if (i,j) is adaption plane 
/COM14/input/ true if (i,k) is adaption plane 
/COM21 /input/ imax for each grid in multi-grid file 
/COM4 /input/ number of points in physical domain, i direction 
/ COM24 / input / IMAX for q file 

/COM21 /input/ if true, describes location to place EXPORT data 
/COM24/ input/ imax for each grid in the multi-grid q file 
/dimension/parameter statement/ currently set at 75 
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INDQ 

INTER 

IQ(NDIM+3) 

IS,IE 

1ST 

ISTEP 

JD 

JEND 

JINVERSE 

JKPLANE 

JM(10) 

JMAX 

JMAXQ 

JMQ(IO) 

JSJE 

JST 

JSTEP 

KD 

KEND 

KINVERSE 

KM(IO) 

KMAX 

KMAXQ 

KMQ(IO) 

KS,KE 

KST 

KSTEP 

LENDADD 

LENDSUB 

LNSING 

LSTADD 

LSTSUB 

MARCH 

MARCHPL 

MAXITS 

MG1 

MG2 

MGCT(IO) 

MGPLS 

MGRID 

MGSTEPS 


/COM5/input/ index for q for adaption variable 
/COM 11 /input/ order of interpolation, 2, 3 or 4 
/COM5/input/ related to INDQ, combines q adaption function 
/COM23/input/ i range of export/import transfer plane 
/COM4/input/ first node along i in adaption domain 
/COM13 /input/ true for stepping in i direction within the plane 

/dimension /parameter statement/ currently set at 75 
/COM4/input/ last node in; adaption domain 
/COM13/INITIAL/ adaption requested in backward; steps 
/COM14/input/ true if (j,k) is adaption plane 
/COM21 /input/ JMAX for each grid in multi-grid file 
/COM4/input/ number of points in physical domain, ;' direction 
/COM24/ input/ JMAX for q file 

/COM24/input/ JMAX for each grid in the multi-grid q file 
/COM23/input/ ; range of export/import transfer plane 
/COM4/input/ first node in ;’ adaption domain 
/COM13/input/ true for stepping in; direction within the plane 

/ dimension/parameter statement/ currently set at 75 
/COM4/input/ last node in k adaption domain 
/COM13/INTTLAL/ adaption requested in backward k steps 
/COM21 /input/ KMAX for each grid in multi-grid file 
/COM4/ input/ number of points in physical domain, k direction 
/COM24/input/ KMAX for q file 

/COM24/input/ KMAX for each grid in the multi-grid q file 
/COM23/input/ k range of export/import transfer plane 
/COM4/input/ first node in k adaption domain 
/COM13/ true for stepping in k direction within the plane 

/COM19/input/ start of add points range 
/COM18/input/ start of delete points range 
/COM15/input/ flag to indicate a common line for each plane 
/COM19/input/ end of add points range 
/COM18/input/ end of delete points range 

/COM14/input/ true to interpolate up to physical boundary 
/COM14/input/ true to interpolate to final plane 
/COM15/INITIAL/ maximum number of iterations for convergence 
/COM9/INITIAL/ number of merging points for NEDGE at i s , 
/COM9/INITIAL/ number of merging points for NEDGE at i tnd 
/COM16/INITIAL/ counter on number of adaptions per grid 
/COM11 /input/ number of planes before full adaption of 
plane parameters 

/COM12/input/ multiple grid number for adaption 

/COM11 /input/ number of lines before full adaption within plane 
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MPLANE /COM23 /input/ plane number for export or import plane 

MSI,MSJ,MSK /COM21 /SIZING/ contains overall maximum size of multiple grids 
MTCH /argument/INITIAL/ flag to indicate export and import has occurred 

MXFG /COM17/MGWALLS/ location of maximum dg/ds 


NADS 

NDIM 

NEDGE 

NFILT 

NFLAG 

NGRID 

NIPTS 

NITGI 

NITGO 

NITQI 

NITQO 

NOMORE 

NOQ 

NOUP 


/COM4/MAIN/ index on number of adaptions 

/dimension/parameter statement/ used for input q, currently set to 5 
/COM9/input/ initiates side-edge boundary override 
/COM11 /input/ number of passes for smoothing data 
/COM15 /INITIAL/ indicates direction sign for normal vectors 
/COM21 /input/ number of grids in multi-grid file 
/ COM4 /INITIAL/ total number of i points in computation domain 
/COM12/INITIAL,WRITOUT.../ unit number for input grid file 
/COM12/INmAL/VVRITOUT.../ unit number for output grid file 
/COM12/INITIAL,WRrrOUT.../ unit number for input q file 
/COM12/INITIAL,WRITOUT.../ unit number for output q file 
/argument/INITIAL/ no more input datasets remain; end program 
/COM11 /input/ false to indicate no q file will be input 
/COM11 /input/ prevents adaption i.e. ( no update 


OK /COM14/SOLUT/ flag to indicate normal termination 

ORTHE(2) /COM16/input/ false removes orthogonality at outer boundaries 

ORTHS(2) /COM1 6 / input/ false removes orthogonality at initial boundaries 


Tl,P2,T>3 / argument / LAGCOF / coefficients of Lagrange polynomials 

PLSING /COM15/input/ flag to define a singular plane 

PRES(IMX) /local/FBAR/ computed pressure 


Q(IDJD,KD,NDIM) /COM1 / input,UPDATE / flow-field variables 

QFUN /COM11 /input/ false to allow geometry function only 

QP(IMX,IMX,NDIM) / COM23 / STOREX / q data of stored plane for grid transfer 


RDSMAX 

RDSMIN 

RE 

REMOVE 

RSWAP 


/COM5/ input/ relative value of As^ 

/COM5/input/ relative value of As mm 
/COM12/input/ not used, PLOT3D variable 
/ COM 18/ input / number of points to remove at outer boundary 
/argument/INITIAL, OUTPUT/ flag for swapping data 


SAVE 

SING 

SMS(IMX) 

SMSS(IMX) 

SN(IMX) 

SNM(IMX) 

SNMK(IMX) 


/COM14/input/ flag to suppress output 
/ argument /NORMPT/ flag to indicate non-existent normal 
/COM3/FBAR/ midpoints of As 
/COM17/WALLS/ s, for geometry calculation 
/COM3/FBAR/ current s array along j line 
/ COM3 /UPDATE/ converged s array along j-1 line 
/COM3/FBAR/ converged s along; line on k-1 plane 
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SP(IMX) 

SPF(IMX) 

SPPL(IMX) 

SS(IMX) 

SUB 

/COM9/TORSION/ s', intersection of torsion vector and ;' line 
/argument/CSPLIN/ cubic spline coefficients 

/COM9 /TORSION/ s‘ , intersection of plane torsion vector and ; line 
/COM3/FBAR/ initial value of s along; line 
/COM18/input/ request to delete points 

TAU 

TAUPL 

TIME 

TN(2) 

TNM(2) 

TRAT(IMX) 

TWOD 

/ local /SOLUT/ x, torsion parameter within plane 
/ local /SOLUT/ y, torsion parameter between planes 
/COM12/input/ not used, PLOT3D variable 

/COMIO /INITIAL/ proportion of b to u for torsion vectors 
/COMIO/TORCOF/ TN modified for boundaries 
/local/FBAR/ computed temperature ratio 
/COM13/input/ indicates 2-D input files 

WDE 

WDES 

WDS 

WDSS 

WEIGHT(IMX) 

WT(IMX) 

WTSUM 

/COM2/WTEDGE/ weighting factor used when NEDGE set 
/COM6/WTEDGE/ retained WDE 
/COM2/WTEDGE/ as WDE, but at initial boundary 
/COM6/WTEDGE/ retained WDS 
/COM2/LINEl,SOLUT/ w, tension parameter 
/COM6/GETWT/ co ( , correction to <o 
/ local /GETB,LINE 1 , WTEDGE / £l/co 

X(ID,JD,KD) 

XJ0MX,3,3) 

XP(IMX,IMX) 

/COM1 /input/ input grid points, i direction 
/COM8 / BLOCK,INTXYZQ/ X data stored in block 
/COM22/STOREX/ stored X of export transfer data 

Y(ID,JD,KD) 

YJ(IMX,3,3) 

YP(IMX,IMX) 

/COM1 /input/ input grid points,; direction 
/COM8/BLOCK,INTXYZQ/ Y data stored in block 
/COM22/STOREX/ stored Y of export transfer data 

Z(ID,JD,KD) 

ZJ(IMX,3,3) 

ZP(IMX,IMX) 

/COM1 /input/ input grid points, k direction 
/ COM8 / BLOCK,INTXYZQ/ Z data stored in block 
/COM22/STOREX/ stored Z of export transfer data 
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3. EXAMPLES 


This section contains many examples to familiarize the user with the adaptive-grid 
process. Each example includes plots of the initial grid and flow-field contours, the input-control- 
parameter file used to adapt the grid, the resultant adapted grid and a discussion on the choice of 
control parameters. In addition, several cases show the improved flow solution obtained from 
the flow solver using the adapted grid as input. The 2-D examples show the effect of within-plane 
parameters that are also needed for 3-D problems, so 3-D users should also study these examples. 

3.1 Two-Dimensional Examples 

The first set of examples are for 2-D problems. Although the code adapts the data as if it 
were a single plane in a 3-D environment, this is imperceptible to the user. When TWOD=.true. 
on the input parameter file, SAGE will expect the data to be in the 2-D format of PLOT3D. The 
index of (1) on the CLAM and CT parameters refers to within-plane variables. Index (2) is not 
needed for these 2-D cases. If TWOD=.t. on the first pass, it will be assumed for subsequent passes 
in the same run. 

Case 1. Flow in a Supersonic Inlet 

Figure 13(a) shows the 101x79 initial grid (assigned to unit 7) for an inlet flow-field 
problem. Since the size of the grid is too large for the parameter statement (defaulted as 
75x75x75) it is first necessary to change the values of ID, JD, KD and IMX (to 101, 101, 1 and 101 
respectively) and recompile SAGE. However, SAGE will warn the user if this is not done. In this 
problem, flow is from left to right and a shock emanates from the upper-wall corner, reflecting 
off the lower wall. In addition, an expansion fan originates from the downstream upper-wall 
corner and interacts with the reflected shock. An interim computed solution (input on unit 8), 
generated by the flow solver using this initial grid, is shown as density contours in Fig. 13(b). 



Figure 13. Flow in a supersonic inlet, (a) Initial grid ; (b) computed density contours. 

The SAGE code is now run to create a new grid that is more adapted to this solution, i.e., 
the grid points are redistributed with respect to a chosen flow-field variable (in this case, density) 
and output to the file assigned to unit 10. SAGE also interpolates the complete flow-field solution 
onto these new grid points and outputs this file to unit 11. The updated files will subsequently be 
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input into the flow solver to produce more accurate solutions. Several examples of grid adaption 
are given for this inlet problem to demonstrate the effect of varying the control parameters. Each 
example uses the grid and solution files shown in Fig. 13 as input. 

Example 1: Single pass, stepping in j direction. Until the user is familiar with the basic 
parameters, the first attempt should be an input-control file with no parameters; i.e., all default 
values are used (with the exception of TWOD=.true.): 


$namel tivod=.t. $ 

On completion of the execution, the message to unit 6 
contains "ADAPTION 1 COMPLETE". The adapted grid 
using the default input parameters is shown in figure 
14(a). The grid has been more evenly spaced and there is 
a slight clustering of points around the shock on the 
lower-wall boundary but the remaining grid lines are 
not adequately adapted to the flow features. This implies 
that the torsion term (whose magnitude is controlled by 
CLAM(l)) is too large and is overriding the local tension 
term, preventing the tension forces from pulling the 
points to the high-gradient regions. In addition, the 
minimum grid spacing is too large. Based on this 
information, a new input-control file was chosen: 

$ndmel twod=.t.,rdsmin=.25,clam(l)=.001,nedge=l $ 



Figure 14. Adapted grids, 
(a) All default parameters. 


Only three control parameters are input: the minimum 

allowable grid spacing, the magnitude of the torsion 
term (an order of magnitude less than the default 
value), and a request for edge control (which is 
frequently used). The remaining parameters retain 
their default values, signifying that the full grid will be 
adapted, density is the adaption variable, and stepping 
is in the j direction. The result of this adaption is 
shown in figure 14(b). It can be seen that, since 
adaption began along the lower surface ( jst=l is the 
default), the redistributed points cluster around the 
point of reflection on this line. However, points do not 
become sufficiently clustered at the corner shock and at 
the start of the expansion wave region on the upper 
surface. This is to be expected, since the adaption on the 
initial line is a 1-D solution and will pick up features 
quite clearly. However, as stepping continues, the 
features on subsequent lines get "dampened" by the 
torsion (i.e., smoothness) control. 



Figure 14 continued, (b) Adaption 
from bottom to top: rdsmin=.25, 
clam(l)=.001, nedge=l. 


Example 2: Adaption with backward /'steps. This example maintains the same parameters as 
example 1, except that the upper surface is chosen as the initial adaption line, i.e., 

$namel twod=.t. / rdsmin=.25,clam(l)=.001 / nedge=l,jst=79,jend=l $ 
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Figure 14(c) shows the result of this adaption. In this case, the upper surface is more clearly 
adapted, and the point of reflection on the lower surface is more spread out. Comparing these 
two examples indicates quite clearly that the starting line has a strong effect on the resultant 
adapted grid, and for some applications it should be chosen carefully. 

Example 3: Adaption in i direction. This example shows the very different grid generated when 
the adaption is performed by stepping in the i direction ( istep=.true .). Based on the experience 
obtained from the first two examples, adapting from right to left will produce a better grid since 
there are no gradients on line ist=l to adapt to. Hence, the input control file was chosen to 
contain: 

$namel twod=.t.,istep=.f.,ist=l01,iend=l,rdsmin=.25,clam(l)=.001,nedge=l $ 

Figure 14(d) shows the result of this adaption. Not surprisingly, the reflected shock region on the 
lower surface is not "captured" by this adaption, and thus this grid is less suitable than those 
shown above. This demonstrates that the choice of stepping direction is important in producing a 
desirable grid. 




Figure 14 continued, (c) Marching from top to bottom; (d) marching in i, from right to left. 

Example 4: Effect of changing control parameters. This example is actually a collection of 
examples that show the effect of varying the control parameters clam(l), ct(l), rdsmax, and 
rdsmin. Figure 15(a) shows ah adapted grid using "baseline" values. These are the same as the 
default values in the code with the exceptions of clam(l)=.001 and adaption proceeds from top to 
bottom (i.e., jst-79, jend=l). 

$namel twod=.true.,clam(l)=.001,jst=79,jend=l,n£dge=l $ 

Figures 15(b) through 15(f) are grids adapted by changing just one of the baseline parameters. 
Figure 15(b) shows the result with clam(l)=.01, and it is immediately obvious why this code 
default value was not chosen as the baseline value for these comparison cases: clam(l) is too large 
to allow any of the flow features to be "captured" by the adaption. Figure 15(c) shows the case for 
clam(l)=.0001, and this smaller value produces a grid with more points clustered around the 
shocks. Figures 15(d) and 15(e) show the effect of the ct(l) parameter: ct(l)=1.0 in Fig. 15(d), 
emphasizing straightness, and ct=.0 in Fig. 15(e), emphasizing orthogonality. The latter shows 
how the orthogonality term dampens the adaption for this case: we are requesting orthogonality 
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(e) 



(f) 


Figure 15. Effect of control parameters, (a) Baseline: clam(l)=.001,ct(l)=.5,rdsmax=2.0,rdsmin=.5; 
(b) clam(l)=M ; (c) clam(l)=.0001 ; (d) ct(l)=1.0; (e) ct(l)=.0; (f) rdsmax=4.0,rdsmin=.25. 
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to the already parallel i lines. Finally, the effect of changing the minimum and maximum 
allowable grid spacing is shown in Fig. 15(f), where rdsmax=4.0 and rdsmin=.25 are used. This 
mesh size control is only a factor of 2 different from that in the baseline example (Fig. 15(a)), but 
significantly densifies the spacing in the shock regions. 

Example 5: Two-directional adaption. Two-directional adaption is created by adapting in one 
direction (e.g., stepping in j) and then adapting in the other direction (i.e., i), using the first 
adapted grid as input. The resultant grid will depend on the order of the stepping direction. Two 
passes (i.e., two sets of input control parameters) were made to produce the adapted grid shown 
in Fig. 16(a). If both these passes are made during the same execution run (not to be 
recommended for a first attempt), then additional output files will be created on units 12 and 13. 
Data will be output for both passes, unless the SAVE parameter is set to false on the first set of 
control parameters. The control file for this example contains: 

$namel twod=.t.,clam(l)=.001,rdsmax=4.0,rdsmin=.25,jst=79,jend=l,nedge=l $ 

$namel istep=.t.,ist=10l / iend=l,clam(l)=.001,rdsmin=.25 / nedge=l $ 



Figure 16. Two-directional adaption, (a) Marching in j followed by marching in i; 

(b) marching in i followed by marching in j. 

The. first pass steps in the j direction and uses the same input control parameters as the 
adapted grid shown in Fig. 15(f). The second pass steps in i, and uses the same parameters as 
example 3. Note that any parameters that remain unchanged for the two passes still need to be 
redefined in the control file, since the code restores all default parameters at the conclusion of 
each pass (with the exception of twod—.true.). The control file used to produce Fig. 16(b) contains 
the same parameters, but the adaption order is reversed. That is, the first pass steps in the i 
direction and the second pass in the j direction. The difference between Figs. 16(a) and 16(b) is 
obvious and shows that the initial grid completely effects the adaption. Although this example 
shows that both adaptive sequences do adapt the grid quite reasonably, it should be noted that it is 
also possible that one order of adaption will produce a completely unacceptable result, while the 
reverse is quite suitable. 

Example 6: Flow solution using adapted grid. Figure 17 demonstrates the improved flow-field 
solution obtained when the flow solver is rerun using an adapted grid as input. The adapted grid 
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shown in Fig. 17(a) was obtained by choosing the 'best' parameters, based on the experience 
already described: 

$namel twod-.true.,jst=79,jend=l,rdsmax=4.0,rdsmin=.25,clani(l)=.0001,nedge=l $ 

This grid and the interpolated flow-field variables were input to the flow solver, producing the 
flow solution (shown as density contours) seen in Fig. 17(b). The improvement in the resolution 
of the incident and reflected shock is obvious when compared to the initial solution shown in 
Fig. 13(b). 



Figure 17. Solution using adapted grid as input, (a) 'Best' adaption, rdsmax=4.0,rdsmin=.25, 
clam(l)=.0001; (b) solution density contours, using 'best' adapted grid. 


Case 2. Hypersonic Blunt-Body Flow 


The second case contains two examples that demonstrate the effect of the CLAM and CT 
parameters. Figures 18(a) and 18(b) show an initial grid (32x32) and corresponding density flow- 
field contours for a hypersonic blunt body. The j direction marches out from the surface to the 
free-stream boundary, and the i direction marches along the body starting at the lowest point. 
This is a simple one-directional problem with the shock shape aligned with the grid. The outer- 
side boundary is curved (when marching in the i direction) and the default values of clam-.Ol 
and ct=.5 will prevent the adapted grid lines from "turning" sufficiently at this edge, giving either 
a false densing of points at the outer side boundary or, possibly, a critical error message. In this 
situation, two remedies are available. CLAM can be decreased, hence de-emphasizing the effect of 
torsion and thereby allowing the tension force to "pull" the nodes toward the shock wave. 
Alternatively, the orthogonality restraint can be removed by setting ct=1.0. Figure 18(c) shows the 
result of setting clam=.001 and retaining all other parameters at their default value. The adapted 
grid obtained by setting ct-1.0 is very similar and is not shown here. The following are the two 
input-control files: 

$namel twod=.t.,istep=.t.,clam(l)=.001$ and $namel twod=.t.,istep=.t.,ct(l)=1.0 $ 

Both adapted grids show the required result, i.e., the clustering of grid points across the shock. 
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Case 3. Blunt-Body Shock Impingement Problem 

The third case introduces the use of the IQ and ORTHE parameters. Figure 19(a) shows the 
input grid of a cowl/lip shock interaction problem. The j direction marches from the body surface 
to the outer free stream and the i direction marches from the lower point of the sphere. The flow- 
field contours of both density (Fig. 19(b)) and Mach number (Fig. 19(c)) are given, since the 
adaption flow-field parameter is a combination of the two. The blunt body shocks, impinging 
shock, and shear layers represent a complex flow that requires adaption in both directions to 
adequately capture all the flow features. Figure 19(d) shows the adapted grid (output on unit 13) 
produced by the following two-pass control file: 

$namel twod=.t.,istep=.t.,ist=91,iend=l,indq=0,iq(l)=2,iq(7)=l,rdsmin-.25, 
ct-.7,clam(l)=.005,nedge=2 $ 

$namel indq=0,iq(l)=2,iq(7)=l,rdsmin=.3,orthe=.f. $ 

The first pass steps in the i direction, starting from the topmost boundary. The parameter 
indq=0 indicates that IQ will be input, and in this case two-thirds of the density function and one- 
third of the Mach function will be combined to become the adaption variable. As explained in 
case 2, the curvature of the outer boundary requires increasing CT. Setting nedge=2 requests that 
only the side-edge boundary at j=l be treated. The second pass introduces the use of the ORTHE 
variable: stepping occurs from the sphere wall to the outer curved boundary, where the default 
would turn the grid to be normal to this outer boundary. ORTHE overrides this and allows the 
adaption to occur naturally. 



Figure 19. Blunt-body shock impingement problem, (a) Initial grid ; (b) initial density contours; 

(c) initial Mach contours; (d) adapted grid. 
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Case 4. Hypersonic Inlet (Zonal Adaption) 

Figures 20(a) and 20(b) show the initial grid and density contours for a hypersonic cowl/lip 
inlet problem. This is a more complex case and requires dividing the adaption domain into two 
zones: the blunt body region and the rectangular inlet region. The upper wall is the j=l line and 
the outgoing channel region (on the right side of the diagram) is the i=l line. 



Figure 20. Hypersonic inlet, zonal adaption, (a) Initial grid; (b) initial density contours. 

Five adaption passes are required to create the final adapted grid shown in Fig. 20(c). The 

input control file consists of 





Figure 20 continued, (c) Adapted grid. 


$namel twod=.t.,rdsmin=1.0 f rdsmax=1.0,save=.f. $ 
$namel jstep-.f.,rdsmin~.2,nedge-l,save=.f. $ 
$namel ist=70,save=.f. $ 

Snamel jst=32,jend=l,iend=71,rdsmin=.25, 
clam(l)=.02,nedge=l,save=.f. $ 

$namel iend=85,jst=19,jend=l,clam(l)=.002, 
nedge=l,mgsteps=5 $ 

The first pass, with equal maximum and 
minimum spacing, spreads out the i grid lines 
evenly; the original grid points were more densely 
distributed in the curved section, leaving fewer 
grid lines in the inlet area. This is a good example 
of using the program to improve an initial grid, 
with no flow-field adaption involved. The second 
pass steps in the i direction and adapts easily 
throughout the entire grid. The third, fourth, and 


fifth passes perform the adaption in the j direction. The very different flow features in the blunt- 
body region (blunt-body shock) compared to the features in the inlet region (Mach stem and 
reflected shocks) indicate dividing the adaption domain into these two zones: i <70 and i>70. 


Only the blunt-body section ( i>70 ) is adapted in pass three, with all default-control parameters. 
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The adaption in the inlet domain starts at jst=32 in order to pick up the Mach stem along the 
lower wall. After stepping through the triple-point region (the intersection of the cowl shock and 
the reflected normal shock), the redistributed points do not spread out sufficiently, so the 
adaption is stopped at line 19 and a fifth pass is started at line 19 with a decreased value of CLAM. 
Since this pass starts internally, it is necessary to set MGSTEPS in order to merge smoothly from 
the already adapted region. Note: save=.f. has been used on the first four passes to prevent 
excessive file generation. 


Case 5. Subsonic Impinging Jet 

This example is to show the use of the SUB parameter. The original grid, with 231x100 
grid points, and the corresponding Mach contours are shown in Figs. 21(a) and (b). It was 
demonstrated that adapting the grid enabled the user to decrease the density of grid points, and 
thus speed up the flow solver. Adapting the grid and reducing the number of points can be 
accomplished in the same execution run. The first step is to ensure that the dimension 
parameters are increased to ID=231, JD=231, KD=1, and IMX=231. (Note that SAGE will inform the 
user of the correct dimensions, if necessary.) To remove every other point from both coordinate 
directions, the input-control file requires two passes for the adaption: 

$namel twod=.t.,sub=l,rdsmax-4.0,rdsmin=.l,clam(l)=.005,jst=100,jend=l,nedge=l,indq=7 $ 
$namel sub=l,istep=.t.,noup=.t. $ 

Both passes use the SUB parameter, reducing the number of points to a 116x51 grid as 
shown in Fig. 20(c). The example shows that the SUB (as with the ADD parameter) can be used in 
conjunction with an adaption (as in the first pass) or simply as a method to reduce the number of 
points in the original grid (as in the second pass). It should be noted that if the two passes had 
been in reverse order, the input value of jst would need to be modified by the user to jst-51. To 
reduce points in a subset of the adaption domain, LSTSUB and LENDSUB should be used, 
however in this example the default values of the entire adaption domain are appropriate. 


(a) 



Figure 21. Subsonic impingement jet. (a) Initial grid. 
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(b) 




(c) 

Figure 21 continued, (b) Initial mach contours; (c) adapted grid, with reduced points. 

Case 6. Axisymmetric Plume Flow 

This final 2-D case is presented to indicate the powerful effect of the adaptive grid process 
on the final flow-field solution. Figure 22(a) shows the initial grid and the reflected computed 
solution (in Mach contours) of an axisymmetric nozzle-plume flow. This initial solution has not 
developed sufficiently to capture the final flow features: the outer shear layer, barrel shock, Mach 
stem, reflected shock, and the triple-point shear layer. Three iterations (through both adaption 
and flow solver) were made to produce the final adapted grid and solution shown in Fig. 22(b). 
The definition of the flow features is greatly improved. Figure 22(c) shows the accuracy of the 
final solution: the lower picture is a shadowgraph of the actual experiment and is almost 
mirrored by the computed solution. 
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Figure 22. Axisymmetric plume flow. 

(a) Initial grid and flow solution (as Mach contours) showing underdeveloped flow features 
(b) final adapted grid and flow solution after several iterations ; 

(c) comparison of computed solution with experimental shadowgraph. 
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Computational domain 



Supersonic flow 


3.2 Three-Dimensional Examples 

The adaption of 3-D grids is more complicated, not only because of the additional torsion 
control, but because of the increased choice of stepping and marching directions. For 2-D 
adaption, four choices of adaption are available: stepping in the i direction, the i j direction, or both 
(in either order). Each of these options will produce a different adapted grid. For 3-D adaption, 
is also necessary to choose a plane-stepping direction and, m theory, to consider up to three 
possible directional passes. In practice, it has been shown that a one-directional pass will 
frequently provide a sufficiently adapted grid, and that a two-directional pass is the maximum 
required. The following case shows the effect of single- and multidirectional passes. 

Case 7: Tandem Fuel Injectors in a Supersonic Combustor (Rockwell Model). 

Figure 23 represents a two-slot, tandem fuel injector arranged behind a backward facing 
6 r step in a supersonic stream that models 

the Rockwell supersonic combustor. The 
fuel injection through the slot nozzle 
creates a complicated three-dimensional 
flow pattern. Since the outflow is 
supersonic, fuel injection normal to the 
main flow produces strong shock waves 
and streamwise separation in the 
vicinity of the slots. In addition, the 
backward facing step locally creates a 
subsonic flow ahead of the slots. In this 
example, grid adaption is limited to the 
region downstream of the backward 
facing step. The grid and initial solutions 
for the fuel injection problem were 
provided by J. Wang. 

Example 1. One adaption pass showing 3-D effect, this example is used to demonstrate how a 
one-directional adaption (i.e., one pass of the program) changes the grid-point distribution, not 
only on the chosen adaption plane, but also on the two cross planes. The initial grid (80x31x61) 
given in Fig 24(a), shows three selected planes from Wang's initial grid, one in each o e 
coordinate directions, where i=40,j=l and k=30. For clarity, these planes are separated (but with 
their original orientations) and shown in Fig. 24(b). Figure 24(c) presents the Mach contours 
obtained from the solution on the initial grid by the flow-field code, corresponding to each plane. 
The adaption is performed by redistributing the points on the first; plane and then marching in ; 
planes (obtained by setting ikplane=.t.). Within each j plane, stepping is in the i direction (since 
istep=.t.) and the adaption is therefore performed along the k lines. The input parameter file used 
to create the adapted grids shown in Fig. 24(d) is: 

$namel ikplane=.t. / istep=.t.,rdsmax=4.0,rdsmin=.2,indq=7,clam(l)=.0002,clam(2)=.0005,nedge=l $ 

As requested, the ; plane has clearly adapted to the Mach number gradients in the k 
direction, but so has the i plane, since this was the stepping direction within planes. However, 
the curvature of the planes themselves has not changed. The k plane shows a different effect. 


Backward 

facing 

step 


Transverse Injectors 


Figure 23. Staged transverse injectors in a supersonic 
combustor, showing computational domain. 
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points have not moved within the plane, but only up and down (again, the k direction), thus 
changing the curvature of the plane. This example clearly shows that although the adaption is in 
only one direction, all planes are affected. It is thus quite possible for one pass to adequately adapt 
the grid for re-input to the flow-field code. 


j=l plane 


j=l plane 


=40 plane 


k=30 plane 




i=40 plane 


k=30 plane 
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Figure 24. Three planes from initial grid, (a) Actual location, i=40,j=l,k=30; 

(b) same planes, separated for clarity; (c) initial Mach contours; (d) single direction adaption. 

Example 2. Two-directional pass. This combustor problem is one that benefits from a two- 
directional pass. By looking at the constant j plane in Fig. 24(d), it can be seen that a second 
adaption to redistribute the i direction would be appropriate: this can be effected by adapting the j 
planes, stepping in k lines. The second set of adaption parameters are 

$namel ikplane=.t.,kstep=.t.,rdsmax=4.0,rdsmin=.2,indq=7,clam(l)=.0002,clam(2)=.0005,nedge=l $ 







Figure 25(a) shows the result of performing this adaption "on top" of the adaption already shown 
in Fig. 24(d); points have now been redistributed in the i direction as well as the k direction. 
However, it should be noted that the i direction can also be adapted by marching in j within 
constant k planes by using: 

$namel ijplane=.t.,jstep=.t.,rdsmax=4.0 / rdsmin=.2,indq=7 r dam(l)=.0002 / dam(2)=.0005,nedge=l $ 

This adapted grid is shown in Fig. 24(b) and shows a very similar redistribution on the j plane, 
even though this is not the adaption plane. 




-30 plane 


(b) 


Figure 25. Two-diredional adaption, using grid shown in Fig. 23(d) as input, 
(a) Adapting ik planes, stepping in k lines; (b) adapting ij planes, stepping in j. 


Example 3. Torsion parameter between planes, A*. Users familiar with 2-D examples know how 
the torsion parameter (A.) affects the adaption within a plane: this parameter controls the 
magnitude of the torsion term and the larger the value of A,, the smoother the resulting grid, but 
with less grid redistribution. This effect is demonstrated in case 1 in the 2-D section. 
Unfortunately, the "base" value of A varies for each problem: the default value in the code is .01 
but the user may have to change this value by as much as one, two or even more orders of 
magnitude. A* is the analogous torsion parameter between planes, and acts in a similar manner 
to A, but restricts the movement of points from plane to plane to maintain smoothness in the 
cross-planes. If A* =0, there is no torsion control between planes, and planes will be adapted as 
independent entities. (In fact, setting A* = 0 is one way of handling periodic planes.) In the code, 
A* is defaulted to .0001 and this example illustrates the effect of this parameter. 

The same combustor-flow problem from example 1 is used with different input control 

parameters. This time, we are adapting k planes and comparing the effect of A on the crossing j 
planes. Figure 26(a) shows the initial distribution of two k planes (the upper and lower surfaces) 
joined by the first j plane. For this particular application, the first k plane would normally not be 
adapted, since the maintenance of its initial spacing is a high priority. However, to illustrate the 


60 



adaption procedure, in this example the first plane is adapted. The input parameter file giving 
the result shown in Fig. 26(b) increases the value of X* to .001: 

$namel ijplane=.t.,jstep=.t.,clam(l)=.001,clam(2)=.001,rdsmax=4.0,rdsmin=.2,nedge=l $ 

This input requests marching in k planes (since ijplane=.t.) and stepping in the j direction within 

the k plane. The only parameter with any control of the j plane is X* ( clam(2 )) and the / plane 
shows some redistribution of points, but insufficient to reflect the flow features. Compare this 
plane to that in Fig. 26(c), where the adaption was performed with the identical input parameters 
except for a decrease in X* to .0001, i.e. 

$namel ijplane=.t. / jstep=.t.,clam(l)=.001,clam(2)=.0001,rdsmax=:4.0,rdsmin=.2,nedge=l $ 

The effect on the smoothing between planes is very evident and indicates the influence of the 
torsion parameter between planes. 



Figure 26. Effect of varying X*. (a) Initial grid showing j=l plane; (b) adaption with clam(2)=.001; 

(c) adaption with clam(2)=.0001. 


Case 8: NASP 3*D Nozzle simulation 

Example 1. Two-directional adaption. Figure 27 shows the geometry of a wind-tunnel 


experimental model for the National AeroSpace 
Nozzle (SERN). The model is inserted into a 
hypersonic test section with cold air injected at 
supersonic speed through the nozzle. 
Superimposed on the ramp section is an 
outline of the computational grid which is 
detailed in Fig. 28(a). This 3-D flow-field grid 
(41x60x90) defines the nozzle and after-body 
region of the model that is used in the 
computational experiments. Shown in the 
figure is the downstream outflow plane at 
i=41, the lower grid at k=31 (part of which is 
the upper surface of the after-body ramp) and 
the symmetry plane at j=l. Additional patched 


Plane (NASP) called the Single Expansion Ramp 



Figure 27. SERN experimental model and the 
computational space in the plume region. 
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Figure 28. NASP nozzle plume flow, (a) initial grid; (b) initial Mach contours; 

(c) adapted grid ; (d) final solution using adapted grid. 

grids are not shown, but there is a grid boundary that must be matched from j— 41 to j—60 along 
the k=31 plane with the grid that is stored in k=l to k=30, and also a matching i=l plane. 
Adaptions were performed from plane i=2 to plane i=41, stepping in both j and k directions 
within the plane/using the initial solution shown as Mach contours in Fig. 28(b). This complex 
case utilizes many of the available input options and the two-pass adaption parameters used to 
produce the adapted grid in Fig. 28(c) were: 

$namel jkplane=.t.,kstep=.t.,rdsmax=2.5,rdsmin=.25,kst=32,kend=83,indq=7 ,nedge=l, 
clam(l)=.l,clam(2)=.01 r mgsteps=4,mgpls=4,march=.t.,ct(l)=.75,save=.f. $ 

$namel jkplane=.t.,jstep=.t. / rdsmax=1.25,rdsmin=.25,kst=35,indq=7,nedge=2,mgl=8, 
clam(l)=.l,clam(2)=.01,mgpls=4,ct(l)=.75 $ 
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The first set of parameters requests the adaption of each i plane (jkplane=.t .), stepping in 
the k direction (kstep=.t.) within the plane. Mach number ( indq=7 ) is the adaption parameter. 
Adaption begins on line k=31 with the merging technique invoked ( mgsteps=4 ) within the plane. 
This will retain the k=31 line to match the existing lower boundary and ensure that the adaption 
parameters are filtered for the next four adaption lines. Both torsion parameters and the 
directional parameter, ct(l) are larger than their default values since a previous test run showed 
that too much adaption (and thus loss of smoothness) occurred using the default values. 
Another parameter used in this example is MARCH. This is invoked from line kend=83: after 
this line, the flow has no gradient features and is mostly numerical "noise", and MARCH simply 
presents a more attractive grid. Since the adaption algorithm normalizes the weighting function, 
"noise" can produce unnecessary redistribution. Due to adjacent grids, it is not appropriate to 
adapt the i=l plane, thus a smooth transition is required between the first plane and subsequent 
planes. This is controlled by mgpls=4: the default start plane of ist=l was not adapted and the next 
three planes were more gradually adapted than they would have been if MGPLS was not 
requested. Example 2 of this case gives a detailed account of this plane-merging process. The 
SAVE parameter was set to .false, to prevent the output datasets from being written. This should 
be done only when more than one adaption pass is made in the same computer run. 

The second adaption set requests adaption of the same i planes, but stepping in the j 
direction within planes. The first adaption point is k=35. This was chosen to retain the first few 
points in the boundary layer. In addition, the parameters nedge=2 and mgl=8 are used to filter 
this boundary layer spacing into the next eight grid points on the line. The value of nedge was 
changed from 1 to 2 for this pass because there was no need to maintain any spacing at the outer 
edge, where no gradients are found: these points are better used within the body of the grid. 
However, the adaption still left an unnecessary number of points in this outer region, so a final 
pass was made: 

$namel jkplane=.t.,jstep=.t.,noup=.t.,sub=l,lstsub=82,add=l,lstadd=51,lendadd=55 $ 

Here, no adaption is performed ( noup=.t .) and four lines are moved by using the SUB and ADD 
parameters. The final adapted grid is given in Fig. 28(c). This was input to the flow solver, 
producing the final solution shown in Fig. 28(d). 

Example 2. Merging from non-adapted planes to adapted planes. The nature of three- 
dimensional problems frequently necessitates the maintenance of boundaries: either for multiple 
grids or to preserve solid geometry walls and even for maintaining boundary layers. For the 
NASP nozzle case described here, the plane at i-1 must match another plane on the adjacent (not 
shown) grid. There is also a dense grid spacing around the nozzle exit for defining the boundary 
layer region that should also be retained. The simplest way to handle these two situations is to 
not adapt the first plane, but to adapt the contiguous planes in a merging manner. This merging 
is needed to create a smooth transition in the cross plane, not in the adaption plane itself (which 
is handled by MGSTEPS). 

Figure 29(a) shows the initial i=l (and thus the identical i=2 plane) for the 3-D grid shown 
in Fig. 28, and Fig. 29(b) contains the corresponding pressure gradients for the i=2 plane. If the 
input adaption parameter ist=2 is used, with no merging process invoked (i.e., the first plane is 
ignored) the resultant adapted i=2 plane is shown in Fig. 29(c). Although this grid is smooth and 
nicely adapted to the flow, the cross plane, is not (compare Fig. 29(c) with Fig. 29(a)). In addition, 
points have been drastically pulled away from the nozzle region. By invoking the MGPLS option, 
this cross plane unevenness can be dampened. The technique is analogous to that used by 
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MGSTEPS: both X* and C* are modified to increase the restraint of movement of points between 
planes. These values will return to their original input values within a certain number of planes, 
as requested by MGPLS. The following is the input parameter file used to create the adapted grid 
( i=2 plane) shown in Fig. 29(d): 

$namel jkplane=.t.,jstep=.t.,rdsmax=1.25,rdsmin=.25,kst=35,indq=6,nedge=2,mgl=8 / 
clam(l)=.l,dam(2)=.01,ct(l)=.75 r mgpls=4 $ 

The difference between Figs. 29(c) and (d) is striking: the boundary layer spacing is now 
maintained and the adaption process is just beginning. Although mgpls=4 in this example, 
setting mgpls=l will also have considerable effect. By setting MGPLS at all, a request is made to 
not adapt the first requested plane (in this case, ist=l, by default) but to use the grid on that plane 
as a control for the subsequent planes. Thus, mgpls=l, ist=l is not the same as setting ist=2. 



Figure 29. Merging planes using mgpls. (a) Initial grid at i=l (and i-2); (b) initial pressure 
contours; (c) adapted grid at i=2 with ist=2; (d) adapted grid at i=2 with ist=l and mgpls=4. 
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Case 9. Aeroassist Flight Experiment (AFE) vehicle 

This example is the hypersonic, non-equilibrium flow around the forebody of the 
aeroassisted flight experiment (AFE) vehicle. In Fig. 30(a) the initial grid configuration (35x23x49) 
is shown around the body in the form of j planes 2 and 22, and the outflow plane at i=35. The 
Mach contours seen in Fig. 30(b) were computed by the flow-field solver of Palmer (1990) using 
the non-adapted grid. This grid was then adapted with respect to these Mach contours, and the 
resulting grid is seen in Fig. 30(c). The redistribution of points within the blunt-body shock region 
is clearly shown. Due to the smooth shape of the shock, adaption was performed only in one 
direction: marching in j planes and stepping in the i direction within planes. The input data set 
used to create the adapted grid in Fig. 30(c) was 

$namel lnsing=l,ikplane=.t.,istep=.t.,indcj=7,jst=2,jend=22,rdsmax=4.0,rdsmin=.l, 
kst-15,orthe=.f.,nedge=l,clam(l)=.001 $ 

This is the first example to use the LNSING parameter. It is set since all j planes emanate 
from a single line at i=l. (This forebody singular line is used in grid mapping.) However, the 
dataset still contains the grid points for line i=l as if it were a separate line in each plane, lnsing-1 
ensures that the first line in the first plane is adapted, and that all i—1 lines on subsequent planes 
are set identical to this adapted line. The parameter orthe was set to .false., which removes the 
orthogonality constraint at the outflow line of z=35. The default value of .true, would have 
created a false turning of the adaption at this location. This adapted grid was then re-input to 
Palmer's flow-solver which showed a considerably sharper blunt-body shock feature (1990). 



Figure 30. Hypersonic forebody flow (AFE). (a) Initial grid around the forebody; 
(b) initial Mach contours; (c) adapted grid. 
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3.3 Multiple-Grid Examples 

Section 2.4.3 describes some of the complexities of multiple grids. The first case here 
reproduces the sty liz ed example shown in Fig. 11 to illustrate the transfer of planes from one grid 
to another. The second case is a generic NASP configuration. (Davies and Deiwert, 1993) 

Case 10. A Simple 3-D Grid 

Example 1. Single plane transfer. Figure 31 is a reproduction of Fig. 11 from section 1. This 
example shows a matching common plane that is stored separately in each grid, but contains 
identical data. We will define the size of each grid as (i,j,k) of (40x15x30). If the grids had been 
separated as single grids, as shown in Fig. 31(b), the only choice would be to start the adaption 
process at the common plane and to march in opposite directions. This would provide identical 
adaption points for the common plane. For some problems this could be quite acceptable, but 
marching in the same direction more often produces the preferred result. 



To adapt multi-zoned grids as shown in Fig. 31(c), a feature is available to transfer data from one 
grid to another. The adaption process in Fig. 31(c) is to: 

(1) adapt grid 1 up to the common plane (j=15 ); 

(2) transfer data from the last plane in grid 1 to the first plane in grid 2; and 

(3) adapt grid 2, using the merging planes option (i.e., mgpls nonzero) for the first plane. 

This sequence is described by the following input control list: 

$namel mgrid=l,ikplane=.t. $ 

$namel export=.t.,mgrid=l,ikplane=.t.,mplane=15 $ 

$namel import=.t.,mgrid=2,ikplane-.t.,mplane=l $ 

$namel mgrid=2,mgpls=2,ikplane=.t. $ 

If only a subset of plane 15, grid 1 was to be transferred to grid 2 (as shown in Fig. 32), then 
the range would be described by entering the non-default values of parameters zs,ic,fcs,fc£ (i.e., 
those relating to the ikplane) for both grids by 

$namel export=.t.,mgrid=l,ikplane=.t.,mplane=15,is=30 $ 

$namel import=.t.,mgrid=2,ikplane=.t.,mplane=l,ie=10 $ 
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In this case, points 30 <i< 40 and all k points from the 15th 
ikplane in grid 1 are transferred to points 1 < i < 10 in the first 
plane of grid 2. (Note that only is=30 on the export list and ie=10 
on the import list are required since the other range values are 
the same as the default values) 

Note for 2-D files . The same concept holds for 2-D files but care 
must be taken since only a line is being transferred. Because only 
one plane exists, ijplane=.t. and mplane=l. Only is,ie r js,je can be 
used and one of these two sets must have equal first and last 
values. Finally, if the export card is the first card in the adaption 
set, include twod=.t. 

Example 2. Multiple Transfers. It is possible to handle more than one transfer within the same 
run of the code. In the example shown in Fig. 32, it is probable that the lighter shaded domain of 
plane 15, grid 1 matches to a third grid, not shown. For this example, the two sets of transfers can 
be entered sequentially: 

$namel exports. t.,mgrid=l,ikplane=.t.,mplane-15,is=30 $ 

$namel import=.t.,mgrid=2,ikplane=.t.,mplane=l,ie=10 $ 

$namel export-.t.,mgrid=l,ikplane=.t.,mplane=15,ie=30 $ 

$namel import=.t.,mgrid=3,.ikplane=.t.,mplane=l $ 

Case 11. Generic NASP configuration. 

This final case is the result of NASP Government Work Package #20, described by Davies 
and Deiwert (1993). A two-part multiple grid defines a generic NASP vehicle configuration: one 
grid (35x59x51) for the nose region and one grid (also 35x59x51) encompassing the remainder of 
the body. Figure 33(a) shows the vehicle surface and three cross planes (jk planes) and Fig. 33(b) 
shows the corresponding Mach contours that were used as the adaption function. The i direction 
is through the vehicle centerline, j steps away from the body and k marches around the body 
from the lower surface. Adaption begins at the nose and marches plane by plane towards the end 
of the body, transferring data at the matching plane {grid 1, jkplane 35 to grid 2, jkplane 1). The 
input control file used to produce the adapted grid shown in Fig. 33(c) is 

$namel mgrid=l,jkplane=.t.,ist=2,kstep=.t.,kst-51,kend=l,jst=15,indq=7,nedge=l, 
rdsmin=.l,rdsmax=3.0,clam(l)=.001,ct(2)=.8,orthe(l)=.f.,orthe(2)=.f. $ 

$namel export=.t.,mgrid=l,jkplane=.t.,mplane=35 $ 

$namel import=.t.,mgrid-2,jkplane=.t.,mplane=l $ 

$namel mgrid=2,jkplane=.t.,kstep=.t.,kst=51,kend=l,indq=7,mgpls=4,nedge=l, 

rdsmin-.l,rdsmax=3.0,clam(l)=.001,ct(2)=.8,orths(2)=.f.,orthe(l)=.f.,orthe(2)=.f. $ 

Some of these adaption parameters need explaining: ist=2 is the starting plane since the nose at 
ist=l is a plane collapsed to a single point; kst=51 and kend-1 produces the best adaption direction 
for kstep (compared to the result using the defaults kst=l and kend=51); jst-15 preserves the 
boundary layer; and ct(2)=.8 (more straightness for the between-planes torsion parameter) 
produces a smoother adaption when looking at the other planar directions. The export and 
import processes are straight forward. For the adaption of grid 2, mgpls=4 and orths(2)=.f. both 
contribute to maintaining continuity across the matching plane. 



Figure 32. Matching subsets. 
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Figure 33 . Generic NASP vehicle configuration, (a) Initial grid; (b) initial Mach contours; and 

(c) adapted grid 
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